diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/CMakeLists.txt b/Detectors/ITSMFT/ITS/postprocessing/studies/CMakeLists.txt index 361ab4db4fb8e..9794b69631d57 100644 --- a/Detectors/ITSMFT/ITS/postprocessing/studies/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/CMakeLists.txt @@ -30,4 +30,6 @@ o2_target_root_dictionary(ITSPostprocessing HEADERS include/ITSStudies/ITSStudiesConfigParam.h include/ITSStudies/TrackCuts.h include/ITSStudies/TrackMethods.h -LINKDEF src/ITSStudiesLinkDef.h) \ No newline at end of file +LINKDEF src/ITSStudiesLinkDef.h) + +add_subdirectory(macros) diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/include/ITSStudies/TrackExtension.h b/Detectors/ITSMFT/ITS/postprocessing/studies/include/ITSStudies/TrackExtension.h index 2567000746559..fd5b93b0f9509 100644 --- a/Detectors/ITSMFT/ITS/postprocessing/studies/include/ITSStudies/TrackExtension.h +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/include/ITSStudies/TrackExtension.h @@ -24,7 +24,7 @@ class MCKinematicsReader; namespace its::study { using mask_t = o2::dataformats::GlobalTrackID::mask_t; -o2::framework::DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, bool useMC, std::shared_ptr kineReader); +o2::framework::DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, std::shared_ptr kineReader); } // namespace its::study } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/macros/CMakeLists.txt b/Detectors/ITSMFT/ITS/postprocessing/studies/macros/CMakeLists.txt new file mode 100644 index 0000000000000..2d78e4077ec53 --- /dev/null +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/macros/CMakeLists.txt @@ -0,0 +1,16 @@ +# 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. + +# o2_add_test_root_macro( +# PostTrackExtension.C +# PUBLIC_LINK_LIBRARIES ROOT::Hist ROOT::RIO ROOT::Core ROOT::Gpad +# LABELS its-study +# COMPILE_ONLY) diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/macros/PostTrackExtension.notest b/Detectors/ITSMFT/ITS/postprocessing/studies/macros/PostTrackExtension.notest new file mode 100644 index 0000000000000..4a7c9c4159a4b --- /dev/null +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/macros/PostTrackExtension.notest @@ -0,0 +1,629 @@ +// 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 "TStyle.h" +#include "TFile.h" +#include "TError.h" +#include "TColor.h" +#include "TCanvas.h" +#include "TH2D.h" +#include "TF1.h" +#include "TEfficiency.h" +#include "TMarker.h" +#include "TLegend.h" +#include "TTree.h" +#include "TLatex.h" + +#include +#include +#include +#endif + +static constexpr std::array bitPatternsBefore{15, 30, 31, 60, 62, 63, 120, 124, 126}; +static constexpr std::array bitPatternsAfter{31, 47, 61, 62, 63, 79, 94, 95, 111, 121, 122, 123, 124, 125, 126, 127}; +inline bool bitsCleared(uint8_t b, uint8_t a) { return (a == b) ? true : (b & ~(a & b)) != 0; } +static constexpr std::array patternColors = { + kRed, // Red + kBlue, // Blue + kGreen, // Green + kMagenta, // Magenta + kCyan, // Cyan + kOrange, // Orange + kViolet, // Violet + kYellow, // Yellow + kPink, // Pink + kAzure, // Azure + kSpring, // Spring Green + kTeal, // Teal + kBlack, // Black + kGray, // Gray + kOrange + 7, // Light Orange + kBlue - 9 // Light Blue +}; + +// Marker styles +static constexpr std::array patternMarkers = { + 20, // Full circle + 21, // Full square + 22, // Full triangle up + 23, // Full triangle down + 24, // Open circle + 25, // Open square + 26, // Open triangle up + 27, // Open cross + 28, // Star + 29, // Plus sign + 30, // Open diamond + 31, // Full diamond + 32, // Cross + 33, // Circle with cross + 34, // X sign + 35 // Double open cross +}; + +enum Labels : unsigned int { + eAll = 0, + eGood, + eFake, + eFakeBefore, + eFakeAfter, + eFakeMix, + eTopGood, + eBotGood, + eMixGood, + eTopFake, + eBotFake, + eMixFake, + eN, +}; +static const std::array names{ + "ALL #frac{ext trks}{all trks}", + "GOOD #frac{good ext trks}{all ext trks}", + "FAKE #frac{fake trks}{all ext trks}", + "FAKE BF #frac{fake bf trks}{fake ext trks}", + "FAKE AF #frac{fake af trks}{fake ext trks}", + "FAKE MIX #frac{fake mix trks}{fake ext trks}", + // Good Top/Bot/Mix + "TOP #frac{good top ext trks}{good ext trks}", + "BOT #frac{good bot ext trks}{good ext trks}", + "MIX #frac{good mix ext trks}{good ext trks}", + // Fake Top/Bot/Mix + "TOP #frac{fake top ext trks}{fake ext trks}", + "BOT #frac{fake bot ext trks}{fake ext trks}", + "MIX #frac{fake mix ext trks}{fake ext trks}", +}; +static const std::array colors{kBlack, kGreen, kRed, kCyan, kYellow, kAzure, + // Good Top/Bot/Mix + kBlue, kOrange, kPink, + // Fake Top/Bot/Mix + kBlue, kOrange, kPink}; +static const std::array markers{20, 21, 22, 23, 27, 28, + // Good Top/Bot/Mix + 29, 33, 39, + // Fake Top/Bot/Mix + 29, 33, 39}; +static const char* const texPtX = "#it{p}_{T} (GeV/#it{c})"; +static const char* const texEff = "Efficiency"; +static const char* const texPtRes = "#sigma(#Delta#it{p}_{T}/#it{p}_{T})"; +static const char* const texDCAxyRes = "#sigma(DCA_{#it{xy}}) (#mum)"; +static const char* const texDCAzRes = "#sigma(DCA_{#it{z}}) (#mum)"; +static const char* const fitOpt{"QWMER"}; + +void setStyle(); +TEfficiency* makeEff(TFile*, const char* num, const char* den); + +template +void style(T* t, Labels lab, TLegend* leg = nullptr) +{ + t->SetMarkerStyle(markers[lab]); + t->SetMarkerColor(colors[lab]); + t->SetLineColor(colors[lab]); + if (leg) { + leg->AddEntry(t, names[lab]); + } +} + +template +void stylePattern(T* t, int i, TLegend* leg = nullptr, const char* name = nullptr) +{ + t->SetMarkerStyle(patternMarkers[i]); + t->SetMarkerColor(patternColors[i]); + t->SetLineColor(patternColors[i]); + if (leg) { + leg->AddEntry(t, name); + } +} + +void PostTrackExtension(const char* fileName = "TrackExtensionStudy.root") +{ + setStyle(); + + std::unique_ptr fIn{TFile::Open(fileName, "READ")}; + if (!fIn || fIn->IsZombie()) { + Error("", "Cannot open file %s", fileName); + return; + } + + { // Purity & Fake-Rate + auto c = new TCanvas("cPFR", "", 800, 600); + auto h = c->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.35, 0.7, 0.7); + auto eff = fIn->Get("eExtension"); + style(eff, eAll, leg); + eff->Draw("same"); + auto effPurity = fIn->Get("eExtensionPurity"); + style(effPurity, eGood, leg); + effPurity->Draw("same"); + auto effFake = fIn->Get("eExtensionFake"); + style(effFake, eFake, leg); + effFake->Draw("same"); + leg->Draw(); + gPad->SetLogx(); + gPad->SetGrid(); + c->SaveAs("trkExt_purity_fake.pdf"); + } + + { // FAKE-Rate composition + auto c = new TCanvas("cFR", "", 800, 600); + auto h = c->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.35, 0.7, 0.7); + auto effFake = fIn->Get("eExtensionFake"); + style(effFake, eFake, leg); + effFake->Draw("same"); + auto effFakeBf = fIn->Get("eExtensionFakeBefore"); + style(effFakeBf, eFakeBefore, leg); + effFakeBf->Draw("same"); + auto effFakeAf = fIn->Get("eExtensionFakeAfter"); + style(effFakeAf, eFakeAfter, leg); + effFakeAf->Draw("same"); + auto effFakeMi = fIn->Get("eExtensionFakeMix"); + style(effFakeMi, eFakeMix, leg); + effFakeMi->Draw("same"); + leg->Draw(); + gPad->SetLogx(); + gPad->SetGrid(); + c->SaveAs("trkExt_fake.pdf"); + } + + { // GOOD Top/Bot/Mix Purity composition + auto c = new TCanvas("cGC", "", 800, 600); + auto h = c->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.35, 0.7, 0.7); + auto effTop = makeEff(fIn.get(), "eExtensionTopPurity", "eExtensionPurity"); + style(effTop, eTopGood, leg); + effTop->Draw("same"); + auto effBot = makeEff(fIn.get(), "eExtensionBotPurity", "eExtensionPurity"); + style(effBot, eBotGood, leg); + effBot->Draw("same"); + auto effMix = makeEff(fIn.get(), "eExtensionMixPurity", "eExtensionPurity"); + style(effMix, eMixGood, leg); + effMix->Draw("same"); + leg->Draw(); + gPad->SetLogx(); + gPad->SetGrid(); + c->SaveAs("trkExt_good_comp.pdf"); + } + + { // FAKE Top/Bot/Mix composition + auto c = new TCanvas("cFC", "", 800, 600); + auto h = c->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.35, 0.7, 0.7); + auto effTop = fIn->Get("eExtensionTopFake"); + style(effTop, eTopFake, leg); + effTop->Draw("same"); + auto effBot = fIn->Get("eExtensionBotFake"); + style(effBot, eBotFake, leg); + effBot->Draw("same"); + auto effMix = fIn->Get("eExtensionMixFake"); + style(effMix, eMixFake, leg); + effMix->Draw("same"); + leg->Draw(); + gPad->SetLogx(); + gPad->SetGrid(); + c->SaveAs("trkExt_fake_comp.pdf"); + } + + { // Good Patterns + auto c = new TCanvas("cPatGood", "", 3 * 800, 3 * 600); + c->Divide(3, 3); + for (int i{0}; i < (int)bitPatternsBefore.size(); ++i) { + auto p = c->cd(i + 1); + auto h = p->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.60, 0.7, 0.88); + leg->SetNColumns(4); + leg->SetHeader(std::format("BEFORE={:07b} GOOD Pattern AFTER/BEFORE", bitPatternsBefore[i]).c_str()); + for (int j{0}; j < (int)bitPatternsAfter.size(); ++j) { + if (bitsCleared(bitPatternsBefore[i], bitPatternsAfter[j])) { + continue; + } + auto eff = fIn->Get(std::format("eExtensionPatternGood_{:07b}_{:07b}", bitPatternsBefore[i], bitPatternsAfter[j]).c_str()); + stylePattern(eff, j, leg, std::format("{:07b}", bitPatternsAfter[j]).c_str()); + eff->Draw("same"); + } + leg->Draw(); + p->SetLogx(); + p->SetGrid(); + } + c->SaveAs("trkExt_good_pattern_comp.pdf"); + } + + { // Fake Patterns + auto c = new TCanvas("cPatFake", "", 3 * 800, 3 * 600); + c->Divide(3, 3); + for (int i{0}; i < (int)bitPatternsBefore.size(); ++i) { + auto p = c->cd(i + 1); + auto h = p->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.60, 0.7, 0.88); + leg->SetNColumns(4); + leg->SetHeader(std::format("BEFORE={:07b} FAKE Pattern AFTER/BEFORE", bitPatternsBefore[i]).c_str()); + for (int j{0}; j < (int)bitPatternsAfter.size(); ++j) { + if (bitsCleared(bitPatternsBefore[i], bitPatternsAfter[j])) { + continue; + } + auto eff = fIn->Get(std::format("eExtensionPatternFake_{:07b}_{:07b}", bitPatternsBefore[i], bitPatternsAfter[j]).c_str()); + stylePattern(eff, j, leg, std::format("{:07b}", bitPatternsAfter[j]).c_str()); + eff->Draw("same"); + } + leg->Draw(); + p->SetLogx(); + p->SetGrid(); + } + c->SaveAs("trkExt_fake_pattern_comp.pdf"); + } + + { // DCA + auto fGaus = new TF1("fGaus", "gaus", -200., 200.); + auto dcaXYVsPtNo = fIn->Get("hDCAxyVsPtResNormal"); + auto dcaXYVsPtYes = fIn->Get("hDCAxyVsPtResExtended"); + auto dcazVsPtNo = fIn->Get("hDCAzVsPtResNormal"); + auto dcazVsPtYes = fIn->Get("hDCAzVsPtResExtended"); + auto bins = dcazVsPtNo->GetXaxis()->GetXbins(); + auto dcaXYResNo = new TH1F("hDcaxyResNo", "NORMAL;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", bins->GetSize() - 1, bins->GetArray()); + auto dcaXYResYes = new TH1F("hDcaxyResYes", "EXTENDED;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", bins->GetSize() - 1, bins->GetArray()); + auto dcaZResNo = new TH1F("hDcazResNo", "NORMAL;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", bins->GetSize() - 1, bins->GetArray()); + auto dcaZResYes = new TH1F("hDcazResYes", "EXTENDED;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", bins->GetSize() - 1, bins->GetArray()); + TH1* proj; + for (int iPt{1}; iPt <= bins->GetSize(); ++iPt) { + auto ptMin = dcaXYResNo->GetXaxis()->GetBinLowEdge(iPt); + if (ptMin < 0.1) { + continue; + } + float minFit = (ptMin < 1.) ? -200. : -75.; + float maxFit = (ptMin < 1.) ? 200. : 75.; + + proj = dcaXYVsPtNo->ProjectionY(Form("hProjDCAxy_no_%d", iPt), iPt, iPt); + proj->Fit("fGaus", fitOpt, "", minFit, maxFit); + dcaXYResNo->SetBinContent(iPt, fGaus->GetParameter(2)); + dcaXYResNo->SetBinError(iPt, fGaus->GetParError(2)); + + proj = dcaXYVsPtYes->ProjectionY(Form("hProjDCAxy_yes_%d", iPt), iPt, iPt); + proj->Fit("fGaus", fitOpt, "", minFit, maxFit); + dcaXYResYes->SetBinContent(iPt, fGaus->GetParameter(2)); + dcaXYResYes->SetBinError(iPt, fGaus->GetParError(2)); + + proj = dcazVsPtNo->ProjectionY(Form("hProjDCAz_no_%d", iPt), iPt, iPt); + proj->Fit("fGaus", fitOpt, "", minFit, maxFit); + dcaZResNo->SetBinContent(iPt, fGaus->GetParameter(2)); + dcaZResNo->SetBinError(iPt, fGaus->GetParError(2)); + + proj = dcazVsPtYes->ProjectionY(Form("hProjDCAz_yes_%d", iPt), iPt, iPt); + proj->Fit("fGaus", fitOpt, "", minFit, maxFit); + dcaZResYes->SetBinContent(iPt, fGaus->GetParameter(2)); + dcaZResYes->SetBinError(iPt, fGaus->GetParError(2)); + } + + dcaXYResNo->SetLineColor(kRed); + dcaXYResNo->SetMarkerColor(kRed); + dcaXYResYes->SetLineColor(kBlue); + dcaXYResYes->SetMarkerColor(kBlue); + dcaZResNo->SetLineColor(kRed); + dcaZResNo->SetMarkerColor(kRed); + dcaZResYes->SetLineColor(kBlue); + dcaZResYes->SetMarkerColor(kBlue); + + auto c = new TCanvas("cDCA", "", 2 * 800, 600); + c->Divide(2, 1); + c->cd(1); + auto h = gPad->DrawFrame(0.1, 1, 10., 500); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texDCAxyRes); + dcaXYResNo->Draw("SAME"); + dcaXYResYes->Draw("SAME"); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + auto leg = new TLegend(0.20, 0.20, 0.40, 0.40); + leg->AddEntry(dcaXYResNo, "Normal"); + leg->AddEntry(dcaXYResYes, "Extended"); + leg->Draw(); + + c->cd(2); + h = gPad->DrawFrame(0.1, 1, 10., 500); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texDCAzRes); + dcaZResNo->Draw("SAME"); + dcaZResYes->Draw("SAME"); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + + c->SaveAs("trkExt_dca.pdf"); + } + + return; + { // Kinematic variables + auto t = fIn->Get("tree"); + auto c = new TCanvas("cKG", "", 800, 600); + c->Divide(3, 2); + { + auto p = c->cd(1); + p->SetGrid(); + auto h = p->DrawFrame(-.6, 0., .6, 9.); + h->GetXaxis()->SetTitle("#frac{Q^{2}}{p_{T,TRK}}-#frac{Q^{2}}{p_{T,MC}}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getQ2Pt()-mcTrk.getQ2Pt()>>hPtNo(100,-.6,.6)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hPtNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.04, 0.04); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-0.55, 8.2, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getQ2Pt()-mcTrk.getQ2Pt()>>hPtYes(100,-.6,.6)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hPtYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.04, 0.04); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-0.55, 7, Form("#mu = %.4f, #sigma = %.4f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(2); + p->SetGrid(); + auto h = p->DrawFrame(-3, 0., 3, 2.); + h->GetXaxis()->SetTitle("Y_{TRK}-Y_{MC}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getY()-mcTrk.getY()>>hYNo(100,-3,3)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hYNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.5, 0.5); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-2, 1.7, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getY()-mcTrk.getY()>>hYYes(100,-3,3)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hYYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.5, 0.5); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-2, 1.5, Form("#mu = %.4f, #sigma = %.4f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(3); + p->SetGrid(); + auto h = p->DrawFrame(-2, 0., 2, 4.2); + h->GetXaxis()->SetTitle("Z_{TRK}-Z_{MC}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getZ()-mcTrk.getZ()>>hZNo(100,-2,2)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hZNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.2, 0.2); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-1.7, 3.8, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getZ()-mcTrk.getZ()>>hZYes(100,-2,2)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hZYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.2, 0.2); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-1.7, 3.5, Form("#mu = %.4f, #sigma = %.4f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(4); + p->SetGrid(); + auto h = p->DrawFrame(-0.02, 0., 0.02, 370.); + h->GetXaxis()->SetTitle("TGL_{TRK}-TGL_{MC}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getTgl()-mcTrk.getTgl()>>hTglNo(100,-0.02,0.02)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hTglNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.003, 0.003); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-0.018, 330, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getTgl()-mcTrk.getTgl()>>hTglYes(100,-0.02,0.02)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hTglYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.003, 0.003); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-0.018, 310, Form("#mu = %.6f, #sigma = %.6f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(5); + p->SetGrid(); + auto h = p->DrawFrame(-0.08, 0., 0.08, 80.); + h->GetXaxis()->SetTitle("SNP_{TRK}-SNP_{MC}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getSnp()-mcTrk.getSnp()>>hSnpNo(100,-0.08,0.08)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hSnpNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.03, 0.03); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-0.07, 72, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getSnp()-mcTrk.getSnp()>>hSnpYes(100,-0.08,0.08)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hSnpYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.03, 0.03); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-0.07, 66, Form("#mu = %.6f, #sigma = %.6f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(6); + auto legend = new TLegend(0.2, 0.2, 0.8, 0.8); + legend->SetTextSize(0.06); + legend->SetLineWidth(3); + legend->SetHeader("GOOD tracks", "C"); + auto mBlue = new TMarker(); + mBlue->SetMarkerColor(kBlue); + mBlue->SetMarkerSize(4); + legend->AddEntry(mBlue, "extended", "p"); + auto mRed = new TMarker(); + mRed->SetMarkerColor(kRed); + mRed->SetMarkerSize(4); + legend->AddEntry(mRed, "normal", "p"); + legend->SetLineColor(kRed); + legend->Draw(); + } + c->SaveAs("trkExt_kinematics.pdf"); + } +} + +void setStyle() +{ + gStyle->Reset("Plain"); + gStyle->SetOptTitle(0); + gStyle->SetOptStat(0); + gStyle->SetPalette(kRainbow); + gStyle->SetCanvasColor(10); + gStyle->SetCanvasBorderMode(0); + gStyle->SetFrameLineWidth(1); + gStyle->SetFrameFillColor(kWhite); + gStyle->SetPadColor(10); + gStyle->SetPadTickX(1); + gStyle->SetPadTickY(1); + gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadLeftMargin(0.15); + gStyle->SetHistLineWidth(1); + gStyle->SetHistLineColor(kRed); + gStyle->SetFuncWidth(2); + gStyle->SetFuncColor(kGreen); + gStyle->SetLineWidth(2); + gStyle->SetLabelSize(0.045, "xyz"); + gStyle->SetLabelOffset(0.01, "y"); + gStyle->SetLabelOffset(0.01, "x"); + gStyle->SetLabelColor(kBlack, "xyz"); + gStyle->SetTitleSize(0.05, "xyz"); + gStyle->SetTitleOffset(1.25, "y"); + gStyle->SetTitleOffset(1.2, "x"); + gStyle->SetTitleFillColor(kWhite); + gStyle->SetTextSizePixels(26); + gStyle->SetTextFont(42); + gStyle->SetTickLength(0.04, "X"); + gStyle->SetTickLength(0.04, "Y"); + gStyle->SetLegendBorderSize(0); + gStyle->SetLegendFillColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetLegendFont(42); +} + +TEfficiency* makeEff(TFile* fIn, const char* num, const char* den) +{ + auto h1 = fIn->Get(num)->GetPassedHistogram(); + auto h2 = fIn->Get(den)->GetPassedHistogram(); + auto e = new TEfficiency(*h1, *h2); + return e; +} diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx b/Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx index 364a354c700b6..465365ffa3d86 100644 --- a/Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx @@ -14,17 +14,24 @@ #include "DataFormatsITS/TrackITS.h" #include "DataFormatsITSMFT/CompCluster.h" #include "DetectorsBase/GRPGeomHelper.h" +#include "DetectorsBase/Propagator.h" #include "Framework/ConfigParamRegistry.h" #include "Framework/Task.h" #include "ITSBase/GeometryTGeo.h" #include "ITSStudies/Helpers.h" #include "ITSStudies/TrackExtension.h" +#include "SimulationDataFormat/MCEventHeader.h" #include "SimulationDataFormat/MCTrack.h" #include "Steer/MCKinematicsReader.h" +#include "ReconstructionDataFormats/Vertex.h" +#include "ReconstructionDataFormats/DCA.h" #include #include "TFile.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TEfficiency.h" namespace o2::its::study { @@ -36,7 +43,9 @@ using o2::steer::MCKinematicsReader; class TrackExtensionStudy : public Task { struct ParticleInfo { - int event; + float eventX; + float eventY; + float eventZ; int pdg; float pt; float eta; @@ -60,24 +69,24 @@ class TrackExtensionStudy : public Task public: TrackExtensionStudy(std::shared_ptr dr, mask_t src, - bool useMC, std::shared_ptr kineReader, std::shared_ptr gr) : mDataRequest(dr), mTracksSrc(src), mKineReader(kineReader), mGGCCDBRequest(gr) { - if (useMC) { - LOGP(info, "Read MCKine reader with {} sources", mKineReader->getNSources()); - } + LOGP(info, "Read MCKine reader with {} sources", mKineReader->getNSources()); } ~TrackExtensionStudy() final = default; void init(InitContext& /*ic*/) final; void run(ProcessingContext& /*pc*/) final; void endOfStream(EndOfStreamContext& /*ec*/) final; + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final; void process(); private: static constexpr std::array mBitPatternsBefore{15, 30, 31, 60, 62, 63, 120, 124, 126}; static constexpr std::array mBitPatternsAfter{31, 47, 61, 62, 63, 79, 94, 95, 111, 121, 122, 123, 124, 125, 126, 127}; + const std::bitset<7> mTopMask{"1110000"}; + const std::bitset<7> mBotMask{"0000111"}; void updateTimeDependentParams(ProcessingContext& pc); std::string mOutFileName = "TrackExtensionStudy.root"; @@ -101,12 +110,30 @@ class TrackExtensionStudy : public Task bool mWithTree{false}; std::unique_ptr mHTrackCounts; - std::unique_ptr mHClustersCounts; std::unique_ptr mHLengthAny, mHLengthGood, mHLengthFake; std::unique_ptr mHChi2Any, mHChi2Good, mHChi2Fake; std::unique_ptr mHPtAny, mHPtGood, mHPtFake; std::unique_ptr mHExtensionAny, mHExtensionGood, mHExtensionFake; - std::unique_ptr mHExtensionPatternsAny, mHExtensionPatternsGood, mHExtensionPatternsFake; + std::unique_ptr mHExtensionPatternsAny, mHExtensionPatternsGood, mHExtensionPatternsFake, mHExtensionPatternsGoodMissed, mHExtensionPatternsGoodEmpty; + std::unique_ptr mEExtensionNum, mEExtensionDen, mEExtensionPurityNum, mEExtensionPurityDen, mEExtensionFakeNum, mEExtensionFakeDen; + std::unique_ptr mEExtensionFakeBeforeNum, mEExtensionFakeAfterNum, mEExtensionFakeMixNum; + std::unique_ptr mEExtensionTopNum, mEExtensionTopPurityNum, mEExtensionTopFakeNum; + std::unique_ptr mEExtensionBotNum, mEExtensionBotPurityNum, mEExtensionBotFakeNum; + std::unique_ptr mEExtensionMixNum, mEExtensionMixPurityNum, mEExtensionMixFakeNum; + std::array, mBitPatternsBefore.size()> mEExtensionPatternGoodNum, mEExtensionPatternFakeNum; + std::array, mBitPatternsAfter.size()>, mBitPatternsBefore.size()> mEExtensionPatternIndGoodNum, mEExtensionPatternIndFakeNum; + // DCA + std::unique_ptr mDCAxyVsPtPionsNormal, mDCAxyVsPtPionsExtended; + std::unique_ptr mDCAzVsPtPionsNormal, mDCAzVsPtPionsExtended; + + template + std::unique_ptr createHistogram(C... n, F... b) + { + auto t = std::make_unique(n..., b...); + mHistograms.push_back(static_cast(t.get())); + return std::move(t); + } + std::vector mHistograms; }; void TrackExtensionStudy::init(InitContext& ic) @@ -114,13 +141,13 @@ void TrackExtensionStudy::init(InitContext& ic) o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); mWithTree = ic.options().get("with-tree"); - constexpr size_t effHistBins = 100; + constexpr size_t effHistBins = 40; constexpr float effPtCutLow = 0.01; constexpr float effPtCutHigh = 10.; auto xbins = helpers::makeLogBinning(effHistBins, effPtCutLow, effPtCutHigh); // Track Counting - mHTrackCounts = std::make_unique("hTrackCounts", "Track Stats", 10, 0, 10); + mHTrackCounts = createHistogram("hTrackCounts", "Track Stats", 10, 0, 10); mHTrackCounts->GetXaxis()->SetBinLabel(1, "Total Tracks"); mHTrackCounts->GetXaxis()->SetBinLabel(2, "Normal ANY Tracks"); mHTrackCounts->GetXaxis()->SetBinLabel(3, "Normal GOOD Tracks"); @@ -132,49 +159,87 @@ void TrackExtensionStudy::init(InitContext& ic) mHTrackCounts->GetXaxis()->SetBinLabel(9, "Extended FAKE AFTER Tracks"); mHTrackCounts->GetXaxis()->SetBinLabel(10, "Extended FAKE BEFORE&AFTER Tracks"); - // Cluster Counting - mHClustersCounts = std::make_unique("hClusterCounts", "Cluster Stats", 5, 0, 5); - mHClustersCounts->GetXaxis()->SetBinLabel(1, "Total Clusters"); - mHClustersCounts->GetXaxis()->SetBinLabel(2, "Tracking"); - mHClustersCounts->GetXaxis()->SetBinLabel(3, "Extension"); - mHClustersCounts->GetXaxis()->SetBinLabel(4, "Good Extension"); - mHClustersCounts->GetXaxis()->SetBinLabel(5, "Fake Extension"); - // Length - mHLengthAny = std::make_unique("hLengthAny", "Extended Tracks Length (ANY);NCluster;Entries", 5, 3, 8); - mHLengthGood = std::make_unique("hLengthGood", "Extended Tracks Length (GOOD);NCluster;Entries", 5, 3, 8); - mHLengthFake = std::make_unique("hLengthFake", "Extended Tracks Length (FAKE);NCluster;Entries", 5, 3, 8); + mHLengthAny = createHistogram("hLengthAny", "Extended Tracks Length (ANY);NCluster;Entries", 5, 3, 8); + mHLengthGood = createHistogram("hLengthGood", "Extended Tracks Length (GOOD);NCluster;Entries", 5, 3, 8); + mHLengthFake = createHistogram("hLengthFake", "Extended Tracks Length (FAKE);NCluster;Entries", 5, 3, 8); // Chi2 - mHChi2Any = std::make_unique("hChi2Any", "Extended Tracks Length (ANY);#chi^{2};Entries", 50, 0, 100); - mHChi2Good = std::make_unique("hChi2Good", "Extended Tracks Length (GOOD);#chi^{2};Entries", 50, 0, 100); - mHChi2Fake = std::make_unique("hChi2Fake", "Extended Tracks Length (FAKE);#chi^{2};Entries", 50, 0, 100); + mHChi2Any = createHistogram("hChi2Any", "Extended Tracks Length (ANY);#chi^{2};Entries", 50, 0, 100); + mHChi2Good = createHistogram("hChi2Good", "Extended Tracks Length (GOOD);#chi^{2};Entries", 50, 0, 100); + mHChi2Fake = createHistogram("hChi2Fake", "Extended Tracks Length (FAKE);#chi^{2};Entries", 50, 0, 100); // Pt - mHPtAny = std::make_unique("hPtAny", "Extended Tracks Length (ANY);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh); - mHPtGood = std::make_unique("hPtGood", "Extended Tracks Length (GOOD);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh); - mHPtFake = std::make_unique("hPtFake", "Extended Tracks Length (FAKE);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh); + mHPtAny = createHistogram("hPtAny", "Extended Tracks Length (ANY);#it{p}_{T};Entries", effHistBins, xbins.data()); + mHPtGood = createHistogram("hPtGood", "Extended Tracks Length (GOOD);#it{p}_{T};Entries", effHistBins, xbins.data()); + mHPtFake = createHistogram("hPtFake", "Extended Tracks Length (FAKE);#it{p}_{T};Entries", effHistBins, xbins.data()); // Length - mHExtensionAny = std::make_unique("hExtensionAny", "Extended Tracks Length (ANY);Extended Layer;Entries", 7, 0, 7); - mHExtensionGood = std::make_unique("hExtensionGood", "Extended Tracks Length (GOOD);Extended Layer;Entries", 7, 0, 7); - mHExtensionFake = std::make_unique("hExtensionFake", "Extended Tracks Length (FAKE);Extended Layer;Entries", 7, 0, 7); + mHExtensionAny = createHistogram("hExtensionAny", "Extended Tracks Length (ANY);Extended Layer;Entries", 7, 0, 7); + mHExtensionGood = createHistogram("hExtensionGood", "Extended Tracks Length (GOOD);Extended Layer;Entries", 7, 0, 7); + mHExtensionFake = createHistogram("hExtensionFake", "Extended Tracks Length (FAKE);Extended Layer;Entries", 7, 0, 7); // Patterns - auto makePatternAxisLabels = [&](TH1* h) { + auto makePatternAxisLabels = [&](TH1* h, bool xBefore = true) { for (int i{1}; i <= h->GetXaxis()->GetNbins(); ++i) { - h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsBefore[i - 1]).c_str()); + if (xBefore) { + h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsBefore[i - 1]).c_str()); + } else { + h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsAfter[i - 1]).c_str()); + } } for (int i{1}; i <= h->GetYaxis()->GetNbins(); ++i) { h->GetYaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsAfter[i - 1]).c_str()); } }; - mHExtensionPatternsAny = std::make_unique("hExtensionPatternsAny", "Extended Tracks Pattern (ANY);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + mHExtensionPatternsAny = createHistogram("hExtensionPatternsAny", "Extended Tracks Pattern (ANY);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); makePatternAxisLabels(mHExtensionPatternsAny.get()); - mHExtensionPatternsGood = std::make_unique("hExtensionPatternsGood", "Extended Tracks Pattern (GOOD);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + mHExtensionPatternsGood = createHistogram("hExtensionPatternsGood", "Extended Tracks Pattern (GOOD);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); makePatternAxisLabels(mHExtensionPatternsGood.get()); - mHExtensionPatternsFake = std::make_unique("hExtensionPatternsFake", "Extended Tracks Pattern (FAKE);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + mHExtensionPatternsFake = createHistogram("hExtensionPatternsFake", "Extended Tracks Pattern (FAKE);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); makePatternAxisLabels(mHExtensionPatternsFake.get()); + mHExtensionPatternsGoodMissed = createHistogram("hExtensionPatternsGoodMissed", "Extended Tracks Pattern (GOOD) Missed Clusters;After;Missed;Entries", mBitPatternsAfter.size(), 0, mBitPatternsAfter.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + makePatternAxisLabels(mHExtensionPatternsGoodMissed.get(), false); + mHExtensionPatternsGoodEmpty = createHistogram("hExtensionPatternsGoodEmpty", "Extended Tracks Pattern (GOOD) Empty Clusters;Before;After;Entries", mBitPatternsAfter.size(), 0, mBitPatternsAfter.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + makePatternAxisLabels(mHExtensionPatternsGoodEmpty.get(), false); + + /// Effiencies + mEExtensionNum = createHistogram("hExtensionNum", "Extension Numerator", effHistBins, xbins.data()); + mEExtensionDen = createHistogram("hExtensionDen", "Extension Dennominator", effHistBins, xbins.data()); + // Purity + mEExtensionPurityNum = createHistogram("hExtensionPurityNum", "Extension Purity Numerator", effHistBins, xbins.data()); + mEExtensionPurityDen = createHistogram("hExtensionPurityDen", "Extension Purity Denominator", effHistBins, xbins.data()); + // Fake + mEExtensionFakeNum = createHistogram("hExtensionFakeNum", "Extension Fake Numerator", effHistBins, xbins.data()); + mEExtensionFakeDen = createHistogram("hExtensionFakeDen", "Extension Fake Denominator", effHistBins, xbins.data()); + mEExtensionFakeBeforeNum = createHistogram("hExtensionFakeBeforeNum", "Extension Fake Before Numerator", effHistBins, xbins.data()); + mEExtensionFakeAfterNum = createHistogram("hExtensionFakeAfterNum", "Extension Fake After Numerator", effHistBins, xbins.data()); + mEExtensionFakeMixNum = createHistogram("hExtensionFakeMixNum", "Extension Fake Mix Numerator", effHistBins, xbins.data()); + // Top + mEExtensionTopNum = createHistogram("hExtensionTopNum", "Extension Top Numerator", effHistBins, xbins.data()); + mEExtensionTopPurityNum = createHistogram("hExtensionTopPurityNum", "Extension Top Purity Numerator", effHistBins, xbins.data()); + mEExtensionTopFakeNum = createHistogram("hExtensionTopFakeNum", "Extension Top Fake Numerator", effHistBins, xbins.data()); + mEExtensionBotNum = createHistogram("hExtensionBotNum", "Extension Bot Numerator", effHistBins, xbins.data()); + mEExtensionBotPurityNum = createHistogram("hExtensionBotPurityNum", "Extension Bot Purity Numerator", effHistBins, xbins.data()); + mEExtensionBotFakeNum = createHistogram("hExtensionBotFakeNum", "Extension Bot Fake Numerator", effHistBins, xbins.data()); + mEExtensionMixNum = createHistogram("hExtensionMixNum", "Extension Mix Numerator", effHistBins, xbins.data()); + mEExtensionMixPurityNum = createHistogram("hExtensionMixPurityNum", "Extension Mix Purity Numerator", effHistBins, xbins.data()); + mEExtensionMixFakeNum = createHistogram("hExtensionMixFakeNum", "Extension Mix Fake Numerator", effHistBins, xbins.data()); + // Patterns + for (int i{0}; i < mBitPatternsBefore.size(); ++i) { + mEExtensionPatternGoodNum[i] = createHistogram(fmt::format("hExtensionPatternGood_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b}", mBitPatternsBefore[i]).c_str(), effHistBins, xbins.data()); + mEExtensionPatternFakeNum[i] = createHistogram(fmt::format("hExtensionPatternFake_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b}", mBitPatternsBefore[i]).c_str(), effHistBins, xbins.data()); + for (int j{0}; j < mBitPatternsAfter.size(); ++j) { + mEExtensionPatternIndGoodNum[i][j] = createHistogram(fmt::format("hExtensionPatternGood_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} -> {:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), effHistBins, xbins.data()); + mEExtensionPatternIndFakeNum[i][j] = createHistogram(fmt::format("hExtensionPatternFake_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} -> {:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), effHistBins, xbins.data()); + } + } + + /// DCA + mDCAxyVsPtPionsNormal = createHistogram("hDCAxyVsPtResNormal", "DCA_{#it{xy}} NORMAL Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500); + mDCAxyVsPtPionsExtended = createHistogram("hDCAxyVsPtResExtended", "DCA_{#it{xy}} EXTENDED Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500); + mDCAzVsPtPionsNormal = createHistogram("hDCAzVsPtResNormal", "DCA_{#it{z}} NORMAL Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500); + mDCAzVsPtPionsExtended = createHistogram("hDCAzVsPtResExtended", "DCA_{#it{z}} EXTENDED Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500); mStream = std::make_unique(mOutFileName.c_str(), "RECREATE"); } @@ -192,8 +257,6 @@ void TrackExtensionStudy::run(ProcessingContext& pc) mClustersMCLCont = recoData.getITSClustersMCLabels(); mInputITSidxs = recoData.getITSTracksClusterRefs(); - mHClustersCounts->SetBinContent(1, mClusters.size()); - LOGP(info, "** Found in {} rofs:\n\t- {} clusters with {} labels\n\t- {} tracks with {} labels", mTracksROFRecords.size(), mClusters.size(), mClustersMCLCont->getIndexedSize(), mTracks.size(), mTracksMCLabels.size()); LOGP(info, "** Found {} sources from kinematic files", mKineReader->getNSources()); @@ -208,10 +271,13 @@ void TrackExtensionStudy::process() for (int iSource{0}; iSource < mKineReader->getNSources(); ++iSource) { mParticleInfo[iSource].resize(mKineReader->getNEvents(iSource)); // events for (int iEvent{0}; iEvent < mKineReader->getNEvents(iSource); ++iEvent) { + const auto& mcEvent = mKineReader->getMCEventHeader(iSource, iEvent); mParticleInfo[iSource][iEvent].resize(mKineReader->getTracks(iSource, iEvent).size()); // tracks for (auto iPart{0}; iPart < mKineReader->getTracks(iEvent).size(); ++iPart) { - auto& part = mKineReader->getTracks(iSource, iEvent)[iPart]; - mParticleInfo[iSource][iEvent][iPart].event = iEvent; + const auto& part = mKineReader->getTracks(iSource, iEvent)[iPart]; + mParticleInfo[iSource][iEvent][iPart].eventX = mcEvent.GetX(); + mParticleInfo[iSource][iEvent][iPart].eventY = mcEvent.GetY(); + mParticleInfo[iSource][iEvent][iPart].eventZ = mcEvent.GetZ(); mParticleInfo[iSource][iEvent][iPart].pdg = part.GetPdgCode(); mParticleInfo[iSource][iEvent][iPart].pt = part.GetPt(); mParticleInfo[iSource][iEvent][iPart].phi = part.GetPhi(); @@ -287,6 +353,8 @@ void TrackExtensionStudy::process() LOGP(info, "\t- Total number of good: {} ({:.2f} %)", good, good * 100. / mTracks.size()); LOGP(info, "\t- Total number of extensions: {} ({:.2f} %)", extended, extended * 100. / mTracks.size()); + o2::dataformats::VertexBase collision; + o2::dataformats::DCA impactParameter; LOGP(info, "** Filling histograms ... "); for (auto iTrack{0}; iTrack < mTracks.size(); ++iTrack) { auto& lab = mTracksMCLabels[iTrack]; @@ -302,6 +370,7 @@ void TrackExtensionStudy::process() continue; } const auto& trk = part.track; + bool isGood = part.isReco && !part.isFake; mHTrackCounts->Fill(0); std::bitset<7> extPattern{0}; @@ -311,44 +380,85 @@ void TrackExtensionStudy::process() } } - if (!extPattern.any()) { - mHTrackCounts->Fill(1); - if (part.isReco || !part.isFake) { - mHTrackCounts->Fill(2); - } else { - mHTrackCounts->Fill(3); + // Tree + while (mWithTree) { + constexpr float refRadius{70.f}; + constexpr float maxSnp{0.9f}; + auto cTrk = trk; + if (!o2::base::Propagator::Instance()->PropagateToXBxByBz(cTrk, refRadius, maxSnp, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrTGeo)) { + break; } - continue; - } - - if (mWithTree) { std::array xyz{(float)part.mcTrack.GetStartVertexCoordinatesX(), (float)part.mcTrack.GetStartVertexCoordinatesY(), (float)part.mcTrack.GetStartVertexCoordinatesZ()}; std::array pxyz{(float)part.mcTrack.GetStartVertexMomentumX(), (float)part.mcTrack.GetStartVertexMomentumY(), (float)part.mcTrack.GetStartVertexMomentumZ()}; auto pdg = O2DatabasePDG::Instance()->GetParticle(part.pdg); if (pdg == nullptr) { - LOGP(fatal, "MC info not available"); + LOGP(error, "MC info not available"); + break; } auto mcTrk = o2::track::TrackPar(xyz, pxyz, TMath::Nint(pdg->Charge() / 3.), true); + if (!mcTrk.rotate(cTrk.getAlpha()) || !o2::base::Propagator::Instance()->PropagateToXBxByBz(mcTrk, refRadius, maxSnp, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrTGeo)) { + break; + } (*mStream) << "tree" - << "trk=" << trk + << "trk=" << cTrk << "mcTrk=" << mcTrk + << "isGood=" << isGood + << "isExtended=" << extPattern.any() << "\n"; + break; + } + + // impact parameter + while (isGood && std::abs(part.pdg) == 211) { + auto trkC = part.track; + collision.setXYZ(part.eventX, part.eventY, part.eventZ); + if (!o2::base::Propagator::Instance()->propagateToDCA(collision, trkC, o2::base::Propagator::Instance()->getNominalBz(), 2.0, o2::base::Propagator::MatCorrType::USEMatCorrTGeo, &impactParameter)) { + break; + } + + auto dcaXY = impactParameter.getY() * 1e4; + auto dcaZ = impactParameter.getZ() * 1e4; + if (!extPattern.any()) { + mDCAxyVsPtPionsNormal->Fill(part.pt, dcaXY); + mDCAzVsPtPionsNormal->Fill(part.pt, dcaZ); + } else { + mDCAxyVsPtPionsExtended->Fill(part.pt, dcaXY); + mDCAzVsPtPionsExtended->Fill(part.pt, dcaZ); + } + break; + } + + mEExtensionDen->Fill(trk.getPt()); + + if (!extPattern.any()) { + mHTrackCounts->Fill(1); + if (part.isReco || !part.isFake) { + mHTrackCounts->Fill(2); + } else { + mHTrackCounts->Fill(3); + } + continue; } mHTrackCounts->Fill(4); mHLengthAny->Fill(trk.getNClusters()); mHChi2Any->Fill(trk.getChi2()); mHPtAny->Fill(trk.getPt()); - if (part.isReco || !part.isFake) { + mEExtensionNum->Fill(trk.getPt()); + mEExtensionPurityDen->Fill(trk.getPt()); + mEExtensionFakeDen->Fill(trk.getPt()); + if (isGood) { mHTrackCounts->Fill(5); mHLengthGood->Fill(trk.getNClusters()); mHChi2Good->Fill(trk.getChi2()); mHPtGood->Fill(trk.getPt()); + mEExtensionPurityNum->Fill(trk.getPt()); } else { mHTrackCounts->Fill(6); mHLengthFake->Fill(trk.getNClusters()); mHChi2Fake->Fill(trk.getChi2()); mHPtFake->Fill(trk.getPt()); + mEExtensionFakeNum->Fill(trk.getPt()); } std::bitset<7> clusPattern{static_cast(trk.getPattern())}; @@ -356,29 +466,28 @@ void TrackExtensionStudy::process() if (extPattern.test(iLayer)) { extPattern.set(iLayer); mHExtensionAny->Fill(iLayer); - if (part.isReco || !part.isFake) { + if (isGood) { mHExtensionGood->Fill(iLayer); } else { mHExtensionFake->Fill(iLayer); } - - if ((part.fakeClusters & (0x1 << iLayer)) == 0) { - mHClustersCounts->Fill(4); - } else { - mHClustersCounts->Fill(5); - } } } - std::bitset<7> oldPattern{clusPattern & ~extPattern}; + std::bitset<7> oldPattern{clusPattern & ~extPattern}, holePattern{clusPattern}; + holePattern.flip(); auto clusN = clusPattern.to_ulong(); auto clusIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), clusN)); auto oldN = oldPattern.to_ulong(); auto oldIdx = std::distance(std::begin(mBitPatternsBefore), std::find(std::begin(mBitPatternsBefore), std::end(mBitPatternsBefore), oldN)); mHExtensionPatternsAny->Fill(oldIdx, clusIdx); - if (part.isReco || !part.isFake) { + if (isGood) { mHExtensionPatternsGood->Fill(oldIdx, clusIdx); + mEExtensionPatternGoodNum[oldIdx]->Fill(trk.getPt()); + mEExtensionPatternIndGoodNum[oldIdx][clusIdx]->Fill(trk.getPt()); } else { mHExtensionPatternsFake->Fill(oldIdx, clusIdx); + mEExtensionPatternFakeNum[oldIdx]->Fill(trk.getPt()); + mEExtensionPatternIndFakeNum[oldIdx][clusIdx]->Fill(trk.getPt()); } // old pattern @@ -392,17 +501,70 @@ void TrackExtensionStudy::process() } } } - if (oldFake && newFake) { mHTrackCounts->Fill(9); + mEExtensionFakeMixNum->Fill(trk.getPt()); } else if (oldFake) { mHTrackCounts->Fill(7); + mEExtensionFakeBeforeNum->Fill(trk.getPt()); } else if (newFake) { mHTrackCounts->Fill(8); + mEExtensionFakeAfterNum->Fill(trk.getPt()); } - mHClustersCounts->SetBinContent(2, mHClustersCounts->GetBinContent(2) + oldPattern.count()); - mHClustersCounts->SetBinContent(3, mHClustersCounts->GetBinContent(3) + extPattern.count()); + // Check if we missed some clusters + if (isGood && holePattern.any()) { + auto missPattern{clusPattern}, emptyPattern{clusPattern}; + for (int iLayer{0}; iLayer < 7; ++iLayer) { + if (!holePattern.test(iLayer)) { + continue; + } + + // Check if there was actually a cluster that we missed + if ((part.clusters & (1 << iLayer)) != 0) { + missPattern.set(iLayer); + } else { + emptyPattern.set(iLayer); + } + } + + if (missPattern != clusPattern) { + auto missN = missPattern.to_ulong(); + auto missIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), missN)); + mHExtensionPatternsGoodMissed->Fill(clusIdx, missIdx); + } + if (emptyPattern != clusPattern) { + auto emptyN = emptyPattern.to_ulong(); + auto emptyIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), emptyN)); + mHExtensionPatternsGoodEmpty->Fill(clusIdx, emptyIdx); + } + } + + // Top/Bot/Mixed Extension + bool isTop = (extPattern & mTopMask).any(); + bool isBot = (extPattern & mBotMask).any(); + if (isTop && isBot) { + mEExtensionMixNum->Fill(trk.getPt()); + if (isGood) { + mEExtensionMixPurityNum->Fill(trk.getPt()); + } else { + mEExtensionMixFakeNum->Fill(trk.getPt()); + } + } else if (isTop) { + mEExtensionTopNum->Fill(trk.getPt()); + if (isGood) { + mEExtensionTopPurityNum->Fill(trk.getPt()); + } else { + mEExtensionTopFakeNum->Fill(trk.getPt()); + } + } else { + mEExtensionBotNum->Fill(trk.getPt()); + if (isGood) { + mEExtensionBotPurityNum->Fill(trk.getPt()); + } else { + mEExtensionBotFakeNum->Fill(trk.getPt()); + } + } } } @@ -421,39 +583,57 @@ void TrackExtensionStudy::endOfStream(EndOfStreamContext& ec) { LOGP(info, "Writing results to {}", mOutFileName); mStream->GetFile()->cd(); + for (const auto h : mHistograms) { + h->Write(); + } - mHTrackCounts->Write(); - mHClustersCounts->Write(); - - mHLengthAny->Write(); - mHLengthGood->Write(); - mHLengthFake->Write(); - - mHChi2Any->Write(); - mHChi2Good->Write(); - mHChi2Fake->Write(); - - mHPtAny->Write(); - mHPtGood->Write(); - mHPtFake->Write(); - - mHExtensionAny->Write(); - mHExtensionGood->Write(); - mHExtensionFake->Write(); - - mHExtensionPatternsAny->Write(); - mHExtensionPatternsGood->Write(); - mHExtensionPatternsFake->Write(); + LOGP(info, "Calculating efficiencies"); + auto makeEff = [](auto num, auto den, const char* name, const char* title) { + auto e = std::make_unique(*num, *den); + e->SetName(name); + e->SetTitle(title); + e->Write(); + }; + makeEff(mEExtensionNum.get(), mEExtensionDen.get(), "eExtension", "Track Extension EXT TRK/ALL"); + makeEff(mEExtensionPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionPurity", "Track Extension Purity GOOD/EXT TRK"); + makeEff(mEExtensionFakeNum.get(), mEExtensionFakeDen.get(), "eExtensionFake", "Track Extension Fake FAKE/EXT TRK"); + makeEff(mEExtensionFakeBeforeNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeBefore", "Track Extension Fake FAKE BEF/FAKE EXT TRK"); + makeEff(mEExtensionFakeAfterNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeAfter", "Track Extension Fake FAKE AFT/FAKE EXT TRK"); + makeEff(mEExtensionFakeMixNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeMix", "Track Extension Fake FAKE MIX/FAKE EXT TRK"); + makeEff(mEExtensionTopNum.get(), mEExtensionDen.get(), "eExtensionTop", "Track Extension Top"); + makeEff(mEExtensionTopPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionTopPurity", "Track Extension Purity GOOD TOP/EXT TRK"); + makeEff(mEExtensionTopFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionTopFake", "Track Extension FAKE TOP/EXT FAKE TRK"); + makeEff(mEExtensionBotNum.get(), mEExtensionDen.get(), "eExtensionBot", "Track Extension Bot"); + makeEff(mEExtensionBotPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionBotPurity", "Track Extension Purity GOOD BOT/EXT TRK"); + makeEff(mEExtensionBotFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionBotFake", "Track Extension FAKE BOT/EXT FAKE TRK"); + makeEff(mEExtensionMixNum.get(), mEExtensionDen.get(), "eExtensionMix", "Track Extension Mix"); + makeEff(mEExtensionMixPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionMixPurity", "Track Extension Purity GOOD MIX/EXT TRK"); + makeEff(mEExtensionMixFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionMixFake", "Track Extension FAKE MIX/EXT FAKE TRK"); + for (int i{0}; i < mBitPatternsBefore.size(); ++i) { + makeEff(mEExtensionPatternGoodNum[i].get(), mEExtensionPurityNum.get(), fmt::format("eExtensionPatternGood_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} GOOD EXT TRK/EXT TRK", mBitPatternsBefore[i]).c_str()); + makeEff(mEExtensionPatternFakeNum[i].get(), mEExtensionFakeNum.get(), fmt::format("eExtensionPatternFake_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} FAKE EXT TRK/EXT TRK", mBitPatternsBefore[i]).c_str()); + for (int j{0}; j < mBitPatternsAfter.size(); ++j) { + makeEff(mEExtensionPatternIndGoodNum[i][j].get(), mEExtensionPatternGoodNum[i].get(), fmt::format("eExtensionPatternGood_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} -> {:07b} GOOD EXT TRK/EXT TRK", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str()); + makeEff(mEExtensionPatternIndFakeNum[i][j].get(), mEExtensionPatternFakeNum[i].get(), fmt::format("eExtensionPatternFake_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} -> {:07b} FAKE EXT TRK/EXT TRK", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str()); + } + } mStream->Close(); } -DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, bool useMC, std::shared_ptr kineReader) +void TrackExtensionStudy::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } +} + +DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, std::shared_ptr kineReader) { std::vector outputs; auto dataRequest = std::make_shared(); - dataRequest->requestTracks(srcTracksMask, useMC); - dataRequest->requestClusters(srcClustersMask, useMC); + dataRequest->requestTracks(srcTracksMask, true); + dataRequest->requestClusters(srcClustersMask, true); auto ggRequest = std::make_shared(false, // orbitResetTime true, // GRPECS=true @@ -468,7 +648,7 @@ DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcCluster "its-study-track-extension", dataRequest->inputs, outputs, - AlgorithmSpec{adaptFromTask(dataRequest, srcTracksMask, useMC, kineReader, ggRequest)}, + AlgorithmSpec{adaptFromTask(dataRequest, srcTracksMask, kineReader, ggRequest)}, Options{{"with-tree", o2::framework::VariantType::Bool, false, {"Produce in addition a tree"}}}}; } diff --git a/Detectors/ITSMFT/ITS/postprocessing/workflow/standalone-postprocessing-workflow.cxx b/Detectors/ITSMFT/ITS/postprocessing/workflow/standalone-postprocessing-workflow.cxx index 02a75def154fc..30fb39c77f235 100644 --- a/Detectors/ITSMFT/ITS/postprocessing/workflow/standalone-postprocessing-workflow.cxx +++ b/Detectors/ITSMFT/ITS/postprocessing/workflow/standalone-postprocessing-workflow.cxx @@ -113,11 +113,14 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.emplace_back(o2::its::study::getAnomalyStudy(srcCls, useMC)); } if (configcontext.options().get("track-extension-study")) { + if (!useMC) { + LOGP(fatal, "Track Extension Study needs MC!"); + } anyStudy = true; srcTrc = GID::getSourcesMask(configcontext.options().get("track-sources")); srcCls = GID::getSourcesMask("ITS"); - o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCls, srcTrc, srcTrc, useMC, srcCls, srcTrc); - specs.emplace_back(o2::its::study::getTrackExtensionStudy(srcTrc, srcCls, useMC, mcKinematicsReader)); + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCls, srcTrc, srcTrc, true, srcCls, srcTrc); + specs.emplace_back(o2::its::study::getTrackExtensionStudy(srcTrc, srcCls, mcKinematicsReader)); } if (configcontext.options().get("efficiency-study")) { anyStudy = true; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index 1c4d3360bc7a3..976d01f1d476b 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -94,11 +94,17 @@ struct TrackingParameters { unsigned long MaxMemory = 12000000000UL; float MaxChi2ClusterAttachment = 60.f; float MaxChi2NDF = 30.f; - bool UseTrackFollower = false; bool FindShortTracks = false; bool PerPrimaryVertexProcessing = false; bool SaveTimeBenchmarks = false; bool DoUPCIteration = false; + /// Cluster attachment + bool UseTrackFollower = false; + bool UseTrackFollowerTop = false; + bool UseTrackFollowerBot = false; + bool UseTrackFollowerMix = false; + float TrackFollowerNSigmaCutZ = 1.f; + float TrackFollowerNSigmaCutPhi = 1.f; }; inline int TrackingParameters::CellMinimumLevel() diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h index 70de43d83d8d2..58483e4aa9f6f 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h @@ -122,7 +122,7 @@ float Tracker::evaluateTask(void (Tracker::*task)(T...), const char* taskName, s { float diff{0.f}; - if (constants::DoTimeBenchmarks) { + if constexpr (constants::DoTimeBenchmarks) { auto start = std::chrono::high_resolution_clock::now(); (this->*task)(std::forward(args)...); auto end = std::chrono::high_resolution_clock::now(); diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h index 207dd5d3d50f5..46499db92d4d5 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h @@ -78,10 +78,10 @@ class TrackerTraits bool isMatLUT() const; // Others - GPUhd() static constexpr int4 getEmptyBinsRect() { return int4{0, 0, 0, 0}; } - const int4 getBinsRect(const Cluster&, int layer, float z1, float z2, float maxdeltaz, float maxdeltaphi); - const int4 getBinsRect(int layer, float phi, float maxdeltaphi, float z, float maxdeltaz); - const int4 getBinsRect(int layer, float phi, float maxdeltaphi, float z1, float z2, float maxdeltaz); + GPUhd() static consteval int4 getEmptyBinsRect() { return int4{0, 0, 0, 0}; } + const int4 getBinsRect(const Cluster&, int layer, float z1, float z2, float maxdeltaz, float maxdeltaphi) const noexcept; + const int4 getBinsRect(int layer, float phi, float maxdeltaphi, float z, float maxdeltaz) const noexcept; + const int4 getBinsRect(int layer, float phi, float maxdeltaphi, float z1, float z2, float maxdeltaz) const noexcept; void SetRecoChain(o2::gpu::GPUChainITS* chain) { mChain = chain; } void setSmoothing(bool v) { mApplySmoothing = v; } bool getSmoothing() const { return mApplySmoothing; } @@ -112,6 +112,12 @@ class TrackerTraits bool mIsGPU = false; }; +inline void TrackerTraits::initialiseTimeFrame(const int iteration) +{ + mTimeFrame->initialise(iteration, mTrkParams[iteration], mTrkParams[iteration].NLayers); + setIsGPU(false); +} + inline float TrackerTraits::getBz() const { return mBz; @@ -122,40 +128,32 @@ inline void TrackerTraits::UpdateTrackingParameters(const std::vectorinitialise(iteration, mTrkParams[iteration], mTrkParams[iteration].NLayers); - setIsGPU(false); -} - -inline const int4 TrackerTraits::getBinsRect(const int layerIndex, float phi, float maxdeltaphi, - float z1, float z2, float maxdeltaz) +inline const int4 TrackerTraits::getBinsRect(const int layerIndex, float phi, float maxdeltaphi, float z1, float z2, float maxdeltaz) const noexcept { const float zRangeMin = o2::gpu::GPUCommonMath::Min(z1, z2) - maxdeltaz; const float phiRangeMin = (maxdeltaphi > constants::math::Pi) ? 0.f : phi - maxdeltaphi; const float zRangeMax = o2::gpu::GPUCommonMath::Max(z1, z2) + maxdeltaz; const float phiRangeMax = (maxdeltaphi > constants::math::Pi) ? constants::math::TwoPi : phi + maxdeltaphi; - if (zRangeMax < -mTrkParams[0].LayerZ[layerIndex + 1] || - zRangeMin > mTrkParams[0].LayerZ[layerIndex + 1] || zRangeMin > zRangeMax) { - + if (zRangeMax < -mTrkParams[0].LayerZ[layerIndex] || + zRangeMin > mTrkParams[0].LayerZ[layerIndex] || zRangeMin > zRangeMax) { return getEmptyBinsRect(); } const IndexTableUtils& utils{mTimeFrame->mIndexTableUtils}; - return int4{o2::gpu::GPUCommonMath::Max(0, utils.getZBinIndex(layerIndex + 1, zRangeMin)), + return int4{o2::gpu::GPUCommonMath::Max(0, utils.getZBinIndex(layerIndex, zRangeMin)), utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMin)), - o2::gpu::GPUCommonMath::Min(mTrkParams[0].ZBins - 1, utils.getZBinIndex(layerIndex + 1, zRangeMax)), // /!\ trkParams can potentially change across iterations + o2::gpu::GPUCommonMath::Min(mTrkParams[0].ZBins - 1, utils.getZBinIndex(layerIndex, zRangeMax)), // /!\ trkParams can potentially change across iterations utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMax))}; } } // namespace its diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h index fe5e52bd6277a..68bfdb51170b5 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h @@ -72,7 +72,9 @@ struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper 0 off + float trackFollowerNSigmaZ = 1.f; // sigma in z-cut for track-following search rectangle + float trackFollowerNSigmaPhi = 1.f; // sigma in phi-cut for track-following search rectangle float cellsPerClusterLimit = -1.f; float trackletsPerClusterLimit = -1.f; int findShortTracks = -1; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h index 7fb5d1421877e..ac0cf51921176 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h @@ -176,7 +176,7 @@ float Vertexer::evaluateTask(void (Vertexer::*task)(T...), const char* taskName, { float diff{0.f}; - if (constants::DoTimeBenchmarks) { + if constexpr (constants::DoTimeBenchmarks) { auto start = std::chrono::high_resolution_clock::now(); (this->*task)(std::forward(args)...); auto end = std::chrono::high_resolution_clock::now(); diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index 9a43402a2e93a..721452bf0361d 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -120,7 +120,7 @@ void Tracker::clustersToTracks(std::function logger, std::f total += evaluateTask(&Tracker::findShortPrimaries, "Short primaries finding", logger); std::stringstream sstream; - if (constants::DoTimeBenchmarks) { + if constexpr (constants::DoTimeBenchmarks) { sstream << std::setw(2) << " - " << "Timeframe " << mTimeFrameCounter++ << " processing completed in: " << total << "ms using " << mTraits->getNThreads() << " threads."; } @@ -200,7 +200,7 @@ void Tracker::clustersToTracksHybrid(std::function logger, // total += evaluateTask(&Tracker::findShortPrimaries, "Hybrid short primaries finding", logger); std::stringstream sstream; - if (constants::DoTimeBenchmarks) { + if constexpr (constants::DoTimeBenchmarks) { sstream << std::setw(2) << " - " << "Timeframe " << mTimeFrameCounter++ << " processing completed in: " << total << "ms using " << mTraits->getNThreads() << " threads."; } @@ -502,8 +502,16 @@ void Tracker::getGlobalConfiguration() if (tc.maxMemory) { params.MaxMemory = tc.maxMemory; } - if (tc.useTrackFollower >= 0) { - params.UseTrackFollower = tc.useTrackFollower; + if (tc.useTrackFollower > 0) { + params.UseTrackFollower = true; + // Bit 0: Allow for mixing of top&bot extension --> implies Bits 1&2 set + // Bit 1: Allow for top extension + // Bit 2: Allow for bot extension + params.UseTrackFollowerMix = ((tc.useTrackFollower & (1 << 0)) != 0); + params.UseTrackFollowerTop = ((tc.useTrackFollower & (1 << 1)) != 0); + params.UseTrackFollowerBot = ((tc.useTrackFollower & (1 << 2)) != 0); + params.TrackFollowerNSigmaCutZ = tc.trackFollowerNSigmaZ; + params.TrackFollowerNSigmaCutPhi = tc.trackFollowerNSigmaPhi; } if (tc.cellsPerClusterLimit >= 0) { params.CellsPerClusterLimit = tc.cellsPerClusterLimit; diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index 4457d4515e0a6..da0abbae9dc1f 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -111,7 +111,7 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROFslice, in const float sqInverseDeltaZ0{1.f / (Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)}; /// protecting from overflows adding the detector resolution const float sigmaZ{o2::gpu::CAMath::Sqrt(Sq(resolution) * Sq(tanLambda) * ((Sq(inverseR0) + sqInverseDeltaZ0) * Sq(meanDeltaR) + 1.f) + Sq(meanDeltaR * tf->getMSangle(iLayer)))}; - const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer, zAtRmin, zAtRmax, + const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer + 1, zAtRmin, zAtRmax, sigmaZ * mTrkParams[iteration].NSigmaCut, tf->getPhiCut(iLayer))}; if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { continue; @@ -679,10 +679,11 @@ void TrackerTraits::extendTracks(const int iteration) for (auto& track : mTimeFrame->getTracks(rof)) { auto backup{track}; bool success{false}; - if (track.getLastClusterLayer() != mTrkParams[iteration].NLayers - 1) { + // the order here biases towards top extension, tracks should probably be fitted separately in the directions and then compared. + if ((mTrkParams[iteration].UseTrackFollowerMix || mTrkParams[iteration].UseTrackFollowerTop) && track.getLastClusterLayer() != mTrkParams[iteration].NLayers - 1) { success = success || trackFollowing(&track, rof, true, iteration); } - if (track.getFirstClusterLayer() != 0) { + if ((mTrkParams[iteration].UseTrackFollowerMix || (mTrkParams[iteration].UseTrackFollowerBot && !success)) && track.getFirstClusterLayer() != 0) { success = success || trackFollowing(&track, rof, false, iteration); } if (success) { @@ -830,8 +831,8 @@ bool TrackerTraits::fitTrack(TrackITSExt& track, int start, int end, int step, f } if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { - float radl = 9.36f; // Radiation length of Si [cm] - float rho = 2.33f; // Density of Si [g/cm^3] + constexpr float radl = 9.36f; // Radiation length of Si [cm] + constexpr float rho = 2.33f; // Density of Si [g/cm^3] if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) { continue; } @@ -855,36 +856,40 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, co auto propInstance = o2::base::Propagator::Instance(); const int step = -1 + outward * 2; const int end = outward ? mTrkParams[iteration].NLayers - 1 : 0; - std::vector hypotheses(1, *track); - for (auto& hypo : hypotheses) { - int iLayer = outward ? track->getLastClusterLayer() : track->getFirstClusterLayer(); + std::vector hypotheses(1, *track); // possibly avoid reallocation + for (size_t iHypo{0}; iHypo < hypotheses.size(); ++iHypo) { + auto hypo{hypotheses[iHypo]}; + int iLayer = static_cast(outward ? hypo.getLastClusterLayer() : hypo.getFirstClusterLayer()); + // per layer we add new hypotheses while (iLayer != end) { - iLayer += step; + iLayer += step; // step through all layers until we reach the end, this allows for skipping on empty layers const float r = mTrkParams[iteration].LayerRadii[iLayer]; + // get an estimate of the trackinf-frame x for the next step float x{-999}; if (!hypo.getXatLabR(r, x, mTimeFrame->getBz(), o2::track::DirAuto) || x <= 0.f) { continue; } - + // estimate hypo's trk parameters at that x auto& hypoParam{outward ? hypo.getParamOut() : hypo.getParamIn()}; if (!propInstance->propagateToX(hypoParam, x, mTimeFrame->getBz(), PropagatorF::MAX_SIN_PHI, PropagatorF::MAX_STEP, mTrkParams[iteration].CorrType)) { continue; } - if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { - float radl = 9.36f; // Radiation length of Si [cm] - float rho = 2.33f; // Density of Si [g/cm^3] + if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { // account for material affects if propagator does not + constexpr float radl = 9.36f; // Radiation length of Si [cm] + constexpr float rho = 2.33f; // Density of Si [g/cm^3] if (!hypoParam.correctForMaterial(mTrkParams[iteration].LayerxX0[iLayer], mTrkParams[iteration].LayerxX0[iLayer] * radl * rho, true)) { continue; } } + + // calculate the search window on this layer const float phi{hypoParam.getPhi()}; const float ePhi{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaSnp2() / hypoParam.getCsp2())}; const float z{hypoParam.getZ()}; const float eZ{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaZ2())}; - const int4 selectedBinsRect{getBinsRect(iLayer, phi, mTrkParams[iteration].NSigmaCut * ePhi, z, mTrkParams[iteration].NSigmaCut * eZ)}; - + const int4 selectedBinsRect{getBinsRect(iLayer, phi, mTrkParams[iteration].TrackFollowerNSigmaCutPhi * ePhi, z, mTrkParams[iteration].TrackFollowerNSigmaCutZ * eZ)}; if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { continue; } @@ -900,9 +905,8 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, co continue; } - TrackITSExt currentHypo{hypo}, newHypo{hypo}; - bool first{true}; - for (int iPhiCount{0}; iPhiCount < phiBinsNum; iPhiCount++) { + // check all clusters in search windows for possible new hypotheses + for (int iPhiCount = 0; iPhiCount < phiBinsNum; iPhiCount++) { int iPhiBin = (selectedBinsRect.y + iPhiCount) % mTrkParams[iteration].PhiBins; const int firstBinIndex{mTimeFrame->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)}; const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1}; @@ -921,7 +925,7 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, co const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer).at(nextCluster.clusterId); - TrackITSExt& tbupdated = first ? hypo : newHypo; + auto tbupdated{hypo}; auto& tbuParams = outward ? tbupdated.getParamOut() : tbupdated.getParamIn(); if (!tbuParams.rotate(trackingHit.alphaTrackingFrame)) { continue; @@ -942,12 +946,7 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, co } tbupdated.setChi2(tbupdated.getChi2() + predChi2); /// This is wrong for outward propagation as the chi2 refers to inward parameters tbupdated.setExternalClusterIndex(iLayer, nextCluster.clusterId, true); - - if (!first) { - hypotheses.emplace_back(tbupdated); - newHypo = currentHypo; - } - first = false; + hypotheses.emplace_back(tbupdated); } } }