diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckTracks.C b/Detectors/ITSMFT/ITS/macros/test/CheckTracks.C index 25a58553a541c..6d9f893225f33 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CheckTracks.C +++ b/Detectors/ITSMFT/ITS/macros/test/CheckTracks.C @@ -71,11 +71,22 @@ void CheckTracks(std::string tracfile = "o2trac_its.root", std::string clusfile Int_t tf = 0, nrec = 0; Int_t lastEventID = -1; Int_t nev = mcTree->GetEntries(); - TH1D* den = new TH1D("den", ";#it{p}_{T} (GeV/#it{c});Den", 100, 0.01, 10); - TH1D* num = new TH1D("num", ";#it{p}_{T} (GeV/#it{c});Num", 100, 0.01, 10); + + Int_t nb = 100; + Double_t xbins[nb + 1], ptcutl = 0.01, ptcuth = 10.; + Double_t a = TMath::Log(ptcuth / ptcutl) / nb; + for (Int_t i = 0; i <= nb; i++) + xbins[i] = ptcutl * TMath::Exp(i * a); + TH1D* num = new TH1D("num", ";#it{p}_{T} (GeV/#it{c});Efficiency (fake-track rate)", nb, xbins); + num->Sumw2(); + TH1D* fak = new TH1D("fak", ";#it{p}_{T} (GeV/#it{c});Fak", nb, xbins); + fak->Sumw2(); + TH1D* den = new TH1D("den", ";#it{p}_{T} (GeV/#it{c});Den", nb, xbins); + den->Sumw2(); + for (Int_t n = 0; n < nev; n++) { std::cout << "\nMC event " << n << '/' << nev << std::endl; - Int_t nGen = 0, nGoo = 0; + Int_t nGen = 0, nGoo = 0, nFak = 0; mcTree->GetEvent(n); Int_t nmc = mcArr->size(); Int_t nmcrefs = mcTrackRefs->size(); @@ -106,7 +117,7 @@ void CheckTracks(std::string tracfile = "o2trac_its.root", std::string clusfile int ok = 0; // Check the availability of clusters - for (int i = 0; i < clusArr->size(); i++) { + for (uint i = 0; i < clusArr->size(); i++) { const Cluster& c = (*clusArr)[i]; auto lab = (clusLabArr->getLabels(i))[0]; if (lab.getEventID() != n) @@ -190,6 +201,9 @@ void CheckTracks(std::string tracfile = "o2trac_its.root", std::string clusfile if (label > 0) { nGoo++; // Good found tracks for the efficiency calculation num->Fill(mcPt); + } else { + nFak++; // Fake-track rate calculation + fak->Fill(mcPt); } } @@ -197,8 +211,11 @@ void CheckTracks(std::string tracfile = "o2trac_its.root", std::string clusfile mcZOut, recZOut, mcPhiOut, recPhiOut, mcThetaOut, recThetaOut, mcPhi, recPhi, mcLam, recLam, mcPt, recPt, ip[0], ip[1], label); } - Float_t eff = (nGen > 0) ? nGoo / Float_t(nGen) : -1.; - std::cout << "Good found tracks: " << nGoo << ", efficiency: " << eff << std::endl; + if (nGen > 0) { + Float_t eff = nGoo / Float_t(nGen); + Float_t rat = nFak / Float_t(nGen); + std::cout << "Good found tracks: " << nGoo << ", efficiency: " << eff << ", fake-track rate: " << rat << std::endl; + } } // "recPt>0" means "found tracks only" @@ -215,9 +232,15 @@ void CheckTracks(std::string tracfile = "o2trac_its.root", std::string clusfile nt->Draw("mcPhiOut-recPhiOut", "recPt>0 && label>0"); new TCanvas; nt->Draw("mcThetaOut-recThetaOut", "recPt>0 && label>0"); - new TCanvas; - num->Divide(den); - num->Draw(); + TCanvas* c1 = new TCanvas; + c1->SetLogx(); + c1->SetGridx(); + c1->SetGridy(); + num->Divide(num, den, 1, 1, "b"); + num->Draw("histe"); + fak->Divide(fak, den, 1, 1, "b"); + fak->SetLineColor(2); + fak->Draw("histesame"); f->Write(); f->Close(); } diff --git a/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CookedTracker.h b/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CookedTracker.h index 7066e134c3dca..237a311bfc81e 100644 --- a/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CookedTracker.h +++ b/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CookedTracker.h @@ -145,7 +145,7 @@ class CookedTracker::Layer void setR(Double_t r) { mR = r; } void unloadClusters(); void selectClusters(std::vector& s, Float_t phi, Float_t dy, Float_t z, Float_t dz); - Int_t findClusterIndex(Double_t z) const; + Int_t findClusterIndex(Float_t z) const; Float_t getR() const { return mR; } const Cluster* getCluster(Int_t i) const { return mClusters[i]; } Float_t getAlphaRef(Int_t i) const { return mAlphaRef[i]; } @@ -161,7 +161,7 @@ class CookedTracker::Layer std::vector mClusters; ///< All clusters std::vector mAlphaRef; ///< alpha of the reference plane std::vector mPhi; ///< cluster phi - std::vector mSectors[kNSectors]; ///< Cluster indices sector-by-sector + std::vector> mSectors[kNSectors]; ///< Cluster indices sector-by-sector }; } // namespace ITS } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/reconstruction/src/CookedTracker.cxx b/Detectors/ITSMFT/ITS/reconstruction/src/CookedTracker.cxx index e80d69906046c..a3fa4caca5f2a 100644 --- a/Detectors/ITSMFT/ITS/reconstruction/src/CookedTracker.cxx +++ b/Detectors/ITSMFT/ITS/reconstruction/src/CookedTracker.cxx @@ -18,6 +18,7 @@ //------------------------------------------------------------------------- #include #include +#include #include #include @@ -45,24 +46,27 @@ using Point3Df = Point3D; // Constants hardcoded for the moment: //************************************************ // seed "windows" in z and phi: makeSeeds -const Double_t kzWin = 0.33; -const Double_t kminPt = 0.05; +const Float_t kzWin = 0.33; +const Float_t kminPt = 0.05; // Maximal accepted impact parameters for the seeds -const Double_t kmaxDCAxy = 3.; -const Double_t kmaxDCAz = 3.; +const Float_t kmaxDCAxy = 3.; +const Float_t kmaxDCAz = 3.; // Layers for the seeding const Int_t kSeedingLayer1 = 6, kSeedingLayer2 = 4, kSeedingLayer3 = 5; // Space point resolution -const Double_t kSigma2 = 0.0005 * 0.0005; +const Float_t kSigma2 = 0.0005 * 0.0005; // Max accepted chi2 -const Double_t kmaxChi2PerCluster = 20.; -const Double_t kmaxChi2PerTrack = 30.; +const Float_t kmaxChi2PerCluster = 20.; +const Float_t kmaxChi2PerTrack = 30.; // Tracking "road" from layer to layer -const Double_t kRoadY = 0.2; -const Double_t kRoadZ = 0.7; +const Float_t kRoadY = 0.2; +const Float_t kRoadZ = 0.3; // Minimal number of attached clusters const Int_t kminNumberOfClusters = 4; +const float kPI = 3.14159f; +const float k2PI = 2 * kPI; + //************************************************ // TODO: //************************************************ @@ -91,10 +95,7 @@ Label CookedTracker::cookLabel(TrackITS& t, Float_t wrong) const // A label<0 indicates that some of the clusters are wrongly assigned. //-------------------------------------------------------------------- Int_t noc = t.getNumberOfClusters(); - std::array lb; - std::array mx; - - int nLabels = 0; + std::map labelOccurence; for (int i = noc; i--;) { Int_t index = t.getClusterIndex(i); @@ -106,28 +107,16 @@ Label CookedTracker::cookLabel(TrackITS& t, Float_t wrong) const if (lab.isEmpty()) break; // all following labels will be empty also // was this label already accounted for ? - bool add = true; - for (int j = nLabels; j--;) { - if (lb[j] == lab) { - add = false; - mx[j]++; // just increment counter - break; - } - } - if (add) { - lb[nLabels] = lab; - mx[nLabels] = 1; - nLabels++; - } + labelOccurence[lab]++; } } Label lab; Int_t maxL = 0; // find most encountered label - for (int i = nLabels; i--;) { - if (mx[i] > maxL) { - maxL = mx[i]; - lab = lb[i]; - } + for (auto[label, count] : labelOccurence) { + if (count <= maxL) + continue; + maxL = count; + lab = label; } if ((1. - Float_t(maxL) / noc) > wrong) { @@ -276,7 +265,7 @@ void CookedTracker::makeSeeds(std::vector& seeds, Int_t first, Int_t l // This is the main pattern recongition function. // Creates seeds out of two clusters and another point. //-------------------------------------------------------------------- - const Double_t zv = getZ(); + const float zv = getZ(); Layer& layer1 = sLayers[kSeedingLayer1]; Layer& layer2 = sLayers[kSeedingLayer2]; @@ -284,6 +273,7 @@ void CookedTracker::makeSeeds(std::vector& seeds, Int_t first, Int_t l const Double_t maxC = TMath::Abs(getBz() * B2C / kminPt); const Double_t kpWin = TMath::ASin(0.5 * maxC * layer1.getR()) - TMath::ASin(0.5 * maxC * layer2.getR()); + const float kpWin100 = kpWin / 100; // Int_t nClusters1 = layer1.getNumberOfClusters(); Int_t nClusters2 = layer2.getNumberOfClusters(); @@ -292,55 +282,54 @@ void CookedTracker::makeSeeds(std::vector& seeds, Int_t first, Int_t l for (Int_t n1 = first; n1 < last; n1++) { const Cluster* c1 = layer1.getCluster(n1); // - // Int_t lab=c1->getLabel(0); + //auto lab = (mClsLabels->getLabels(c1-mFirstCluster))[0]; // - Double_t z1 = c1->getZ(); auto xyz1 = c1->getXYZGloRot(*mGeom); - Double_t r1 = xyz1.rho(), phi1 = layer1.getClusterPhi(n1); + auto z1 = xyz1.Z(); + auto r1 = xyz1.rho(); + + auto phi1 = layer1.getClusterPhi(n1); + auto tgl = std::abs((z1 - zv) / r1); - Double_t zr2 = zv + layer2.getR() / r1 * (z1 - zv); - Int_t start2 = layer2.findClusterIndex(zr2 - kzWin); + auto zr2 = zv + layer2.getR() / r1 * (z1 - zv); + auto phir2 = phi1; + auto dz2 = kzWin * (1 + 2 * tgl); - for (Int_t n2 = start2; n2 < nClusters2; n2++) { + std::vector selected2; + float dy2 = kpWin * layer2.getR(); + layer2.selectClusters(selected2, phir2, dy2, zr2, dz2); + for (auto n2 : selected2) { const Cluster* c2 = layer2.getCluster(n2); // - // if (c2->getLabel(0)!=lab) continue; + //if ((mClsLabels->getLabels(c2-mFirstCluster))[0] != lab) continue; // - Double_t z2 = c2->getZ(); - if (z2 > (zr2 + kzWin)) - break; // check in Z - - Double_t phi2 = layer2.getClusterPhi(n2); - if (TMath::Abs(phi2 - phi1) > kpWin) - continue; // check in Phi - auto xyz2 = c2->getXYZGloRot(*mGeom); - Double_t r2 = xyz2.rho(); - Double_t crv = f1(xyz1.X(), xyz1.Y(), xyz2.X(), xyz2.Y(), getX(), getY()); + auto z2 = xyz2.Z(); + auto r2 = xyz2.rho(); + + Float_t hcrv = 0.5 * f1(xyz1.X(), xyz1.Y(), xyz2.X(), xyz2.Y(), getX(), getY()); - Double_t zr3 = z1 + (layer3.getR() - r1) / (r2 - r1) * (z2 - z1); - Double_t dz = kzWin / 2; + auto zr3 = z1 + (layer3.getR() - r1) / (r2 - r1) * (z2 - z1); + auto phir3 = phi1 + hcrv * (layer3.getR() - r1); + auto dz3 = 0.5f * dz2; - Int_t start3 = layer3.findClusterIndex(zr3 - dz); - for (Int_t n3 = start3; n3 < nClusters3; n3++) { + std::vector selected3; + float dy3 = kpWin100 * layer3.getR(); + layer3.selectClusters(selected3, phir3, dy3, zr3, dz3); + for (auto n3 : selected3) { const Cluster* c3 = layer3.getCluster(n3); // - // if (c3->getLabel(0)!=lab) continue; + //if ((mClsLabels->getLabels(c3-mFirstCluster))[0] != lab) continue; // - Double_t z3 = c3->getZ(); - if (z3 > (zr3 + dz)) - break; // check in Z - - Double_t r3 = c3->getX(); - Double_t phir3 = phi1 + 0.5 * crv * (r3 - r1); - Double_t phi3 = layer3.getClusterPhi(n3); - if (TMath::Abs(phir3 - phi3) > kpWin / 100) - continue; // check in Phi + auto xyz3 = c3->getXYZGloRot(*mGeom); + auto z3 = xyz3.Z(); + auto r3 = xyz3.rho(); - Point3Df txyz2 = c2->getXYZ(); // tracking coordinates - // txyz2.SetX(layer2.getXRef(n2)); // The clusters are already in the tracking frame + zr3 = z1 + (r3 - r1) / (r2 - r1) * (z2 - z1); + if (std::abs(z3 - zr3) > 0.2 * dz3) + continue; - auto xyz3 = c3->getXYZGloRot(*mGeom); + const Point3Df& txyz2 = c2->getXYZ(); // tracking coordinates TrackITS seed = cookSeed(xyz1, xyz3, txyz2, layer2.getR(), layer3.getR(), layer2.getAlphaRef(n2), getBz()); @@ -394,26 +383,22 @@ void CookedTracker::trackSeeds(std::vector& seeds) } for (auto& track : seeds) { - Double_t x = track.getX(); - Double_t y = track.getY(); - Double_t phi = track.getAlpha() + TMath::ATan2(y, x); - const Float_t pi2 = 2. * TMath::Pi(); - if (phi < 0.) - phi += pi2; - else if (phi >= pi2) - phi -= pi2; - - Double_t z = track.getZ(); - Double_t crv = track.getCurvature(getBz()); - Double_t tgl = track.getTgl(); - Double_t r1 = sLayers[kSeedingLayer2].getR(); + auto x = track.getX(); + auto y = track.getY(); + Float_t phi = track.getAlpha() + TMath::ATan2(y, x); + BringTo02Pi(phi); + + auto z = track.getZ(); + auto crv = track.getCurvature(getBz()); + auto tgl = track.getTgl(); + Float_t r1 = sLayers[kSeedingLayer2].getR(); for (Int_t l = kSeedingLayer2 - 1; l >= 0; l--) { - Double_t r2 = sLayers[l].getR(); + Float_t r2 = sLayers[l].getR(); phi += 0.5 * crv * (r2 - r1); z += tgl / (0.5 * crv) * (TMath::ASin(0.5 * crv * r2) - TMath::ASin(0.5 * crv * r1)); selec[l].clear(); - sLayers[l].selectClusters(selec[l], phi, kRoadY, z, kRoadZ); + sLayers[l].selectClusters(selec[l], phi, kRoadY, z, kRoadZ * (1 + 2 * std::abs(tgl))); r1 = r2; } @@ -528,7 +513,7 @@ void CookedTracker::process(const std::vector& clusters, std::vectorgetZ() < c2->getZ()); }); Double_t r = 0.; - const Float_t pi2 = 2. * TMath::Pi(); Int_t m = mClusters.size(); for (Int_t i = 0; i < m; i++) { const Cluster* c = mClusters[i]; @@ -717,8 +701,8 @@ void CookedTracker::Layer::init() Float_t phi = xyz.Phi(); BringTo02Pi(phi); mPhi.push_back(phi); - Int_t s = phi * kNSectors / pi2; - mSectors[s].push_back(i); + Int_t s = phi * kNSectors / k2PI; + mSectors[s].emplace_back(i, c->getZ()); } if (m) @@ -746,13 +730,13 @@ Bool_t CookedTracker::Layer::insertCluster(const Cluster* c) return kTRUE; } -Int_t CookedTracker::Layer::findClusterIndex(Double_t z) const +Int_t CookedTracker::Layer::findClusterIndex(Float_t z) const { //-------------------------------------------------------------------- // This function returns the index of the first cluster with its fZ >= "z". //-------------------------------------------------------------------- auto found = std::upper_bound(std::begin(mClusters), std::end(mClusters), z, - [](Double_t zc, const Cluster* c) { return (zc < c->getZ()); }); + [](Float_t zc, const Cluster* c) { return (zc < c->getZ()); }); return found - std::begin(mClusters); } @@ -764,37 +748,31 @@ void CookedTracker::Layer::selectClusters(std::vector& selec, Float_t phi Float_t zMin = z - dz; Float_t zMax = z + dz; - const Float_t pi2 = 2. * TMath::Pi(); + BringTo02Pi(phi); + Float_t dphi = dy / mR; - Float_t phiMin = phi - dphi; - Float_t phiMax = phi + dphi; - Float_t phiRange[2]{ phiMin, phiMax }; + int smin = (phi - dphi) / k2PI * kNSectors; + int ds = (phi + dphi) / k2PI * kNSectors - smin + 1; - Int_t n = 0; - Int_t sector = -1; - for (auto phiM : phiRange) { - Int_t s = phiM * kNSectors / pi2; - if (s < 0) - s += kNSectors; - else if (s >= kNSectors) - s -= kNSectors; + smin = (smin + kNSectors) % kNSectors; - if (s == sector) - break; - sector = s; + for (int is = 0; is < ds; is++) { + Int_t s = (smin + is) % kNSectors; - auto cmp = [this](Double_t zc, Int_t ic) { return (zc < mClusters[ic]->getZ()); }; + auto cmp = [](Float_t zc, std::pair p) { return (zc < p.second); }; auto imin = std::upper_bound(std::begin(mSectors[s]), std::end(mSectors[s]), zMin, cmp); auto imax = std::upper_bound(imin, std::end(mSectors[s]), zMax, cmp); for (; imin != imax; imin++) { - Int_t i = *imin; - Float_t cphi = mPhi[i]; - if (cphi <= phiMin) - continue; - if (cphi > phiMax) - continue; - + auto[i, zz] = *imin; + auto cdphi = std::abs(mPhi[i] - phi); + if (cdphi > dphi) { + if (cdphi > kPI) { + cdphi = k2PI - cdphi; + } + if (cdphi > dphi) + continue; // check in Phi + } selec.push_back(i); } } diff --git a/Detectors/ITSMFT/ITS/workflow/src/CookedTrackerSpec.cxx b/Detectors/ITSMFT/ITS/workflow/src/CookedTrackerSpec.cxx index 2a06835b22174..f532f5fe1d05c 100644 --- a/Detectors/ITSMFT/ITS/workflow/src/CookedTrackerSpec.cxx +++ b/Detectors/ITSMFT/ITS/workflow/src/CookedTrackerSpec.cxx @@ -37,6 +37,8 @@ namespace ITS void CookedTrackerDPL::init(InitContext& ic) { + auto nthreads = ic.options().get("nthreads"); + mTracker.setNumberOfThreads(nthreads); auto filename = ic.options().get("grp-file"); const auto grp = o2::parameters::GRPObject::loadFrom(filename.c_str()); if (grp) { @@ -110,6 +112,7 @@ DataProcessorSpec getCookedTrackerSpec() AlgorithmSpec{ adaptFromTask() }, Options{ { "grp-file", VariantType::String, "o2sim_grp.root", { "Name of the output file" } }, + { "nthreads", VariantType::Int, 1, { "Number of threads" } }, } }; }