AliceO2Group · davidrohr · Sep 18, 2021 · Jul 21, 2021 · Jul 29, 2021 · Aug 10, 2021
@@ -95,9 +95,10 @@ class TrackITS : public o2::track::TrackParCov
  const o2::track::TrackParCov& getParamOut() const { return mParamOut; }

  void setPattern(uint32_t p) { mPattern = p; }
-  int getPattern() const { return mPattern; }
+  uint32_t getPattern() const { return mPattern; }
  bool hasHitOnLayer(int i) { return mPattern & (0x1 << i); }
  bool isFakeOnLayer(int i) { return !(mPattern & (0x1 << (16 + i))); }
+  int getNFakeClusters();

  void setNextROFbit(bool toggle = true) { setUserField((getUserField() & ~kNextROF) | (-toggle & kNextROF)); }
  bool hasHitInNextROF() const { return getUserField() & kNextROF; }
@@ -132,6 +133,8 @@ class TrackITSExt : public TrackITS
    setNumberOfClusters(ncl);
  }

+  GPUdDefault() TrackITSExt(const TrackITSExt& t) = default;
+
  void setClusterIndex(int l, int i)
  {
    int ncl = getNumberOfClusters();

@@ -106,3 +106,14 @@ Bool_t TrackITS::isBetter(const TrackITS& best, Float_t maxChi2) const
  }
  return kFALSE;
 }
+
+int TrackITS::getNFakeClusters()
+{
+  int nFake{0};
+  for (int iCl{0}; iCl < getNClusters(); ++iCl) {
+    if (hasHitOnLayer(iCl) && isFakeOnLayer(iCl)) {
+      ++nFake;
+    }
+  }
+  return nFake;
+}
@@ -12,6 +12,9 @@
 #include <TClonesArray.h>
 #include <TH2F.h>
 #include <TCanvas.h>
+#include <THStack.h>
+#include <TLegend.h>
+#include <TPad.h>

 #include "ITSBase/GeometryTGeo.h"
 #include "SimulationDataFormat/TrackReference.h"
@@ -43,7 +46,7 @@ struct ParticleInfo {

 #pragma link C++ class ParticleInfo + ;

-void CheckTracksCA(std::string tracfile = "o2trac_its.root", std::string clusfile = "o2clus_its.root", std::string kinefile = "o2sim_Kine.root")
+void CheckTracksCA(bool doFakeClStud = true, std::string tracfile = "o2trac_its.root", std::string clusfile = "o2clus_its.root", std::string kinefile = "o2sim_Kine.root")
 {

  using namespace o2::itsmft;
@@ -271,4 +274,93 @@ void CheckTracksCA(std::string tracfile = "o2trac_its.root", std::string clusfil
  clone->Write("clones");
  file.Close();
  std::cout << " done." << std::endl;
+
+  //////////////////////
+  // Fake clusters study
+  if (doFakeClStud) {
+    std::vector<TH1I*> histLength, histLength1Fake, histLengthNoCl, histLength1FakeNoCl;
+    std::vector<THStack*> stackLength, stackLength1Fake;
+    std::vector<TLegend*> legends, legends1Fake;
+    histLength.resize(4);
+    histLength1Fake.resize(4);
+    histLengthNoCl.resize(4);
+    histLength1FakeNoCl.resize(4);
+    stackLength.resize(4);
+    stackLength1Fake.resize(4);
+    legends.resize(4);
+    legends1Fake.resize(4);
+
+    for (int iH{4}; iH < 8; ++iH) {
+      histLength[iH - 4] = new TH1I(Form("trk_len_%d", iH), "#exists cluster", 7, -.5, 6.5);
+      histLength[iH - 4]->SetFillColor(kBlue);
+      histLength[iH - 4]->SetLineColor(kBlue);
+      histLength[iH - 4]->SetFillStyle(3352);
+      histLengthNoCl[iH - 4] = new TH1I(Form("trk_len_%d_nocl", iH), "#slash{#exists} cluster", 7, -.5, 6.5);
+      histLengthNoCl[iH - 4]->SetFillColor(kRed);
+      histLengthNoCl[iH - 4]->SetLineColor(kRed);
+      histLengthNoCl[iH - 4]->SetFillStyle(3352);
+      stackLength[iH - 4] = new THStack(Form("stack_trk_len_%d", iH), Form("trk_len=%d", iH));
+      stackLength[iH - 4]->Add(histLength[iH - 4]);
+      stackLength[iH - 4]->Add(histLengthNoCl[iH - 4]);
+    }
+    for (int iH{4}; iH < 8; ++iH) {
+      histLength1Fake[iH - 4] = new TH1I(Form("trk_len_%d_1f", iH), "#exists cluster", 7, -.5, 6.5);
+      histLength1Fake[iH - 4]->SetFillColor(kBlue);
+      histLength1Fake[iH - 4]->SetLineColor(kBlue);
+      histLength1Fake[iH - 4]->SetFillStyle(3352);
+      histLength1FakeNoCl[iH - 4] = new TH1I(Form("trk_len_%d_1f_nocl", iH), "#slash{#exists} cluster", 7, -.5, 6.5);
+      histLength1FakeNoCl[iH - 4]->SetFillColor(kRed);
+      histLength1FakeNoCl[iH - 4]->SetLineColor(kRed);
+      histLength1FakeNoCl[iH - 4]->SetFillStyle(3352);
+      stackLength1Fake[iH - 4] = new THStack(Form("stack_trk_len_%d_1f", iH), Form("trk_len=%d, 1 Fake", iH));
+      stackLength1Fake[iH - 4]->Add(histLength1Fake[iH - 4]);
+      stackLength1Fake[iH - 4]->Add(histLength1FakeNoCl[iH - 4]);
+    }
+
+    for (auto& event : info) {
+      for (auto& part : event) {
+        int nCl{0};
+        for (unsigned int bit{0}; bit < sizeof(pInfo.clusters) * 8; ++bit) {
+          nCl += bool(part.clusters & (1 << bit));
+        }
+        if (nCl < 3) {
+          continue;
+        }
+
+        auto& track = part.track;
+        auto len = track.getNClusters();
+        for (int iLayer{0}; iLayer < 7; ++iLayer) {
+          if (track.hasHitOnLayer(iLayer)) {
+            if (track.isFakeOnLayer(iLayer)) {       // Reco track has fake cluster
+              if (part.clusters & (0x1 << iLayer)) { // Correct cluster exists
+                histLength[len - 4]->Fill(iLayer);
+                if (track.getNFakeClusters() == 1) {
+                  histLength1Fake[len - 4]->Fill(iLayer);
+                }
+              } else {
+                histLengthNoCl[len - 4]->Fill(iLayer);
+                if (track.getNFakeClusters() == 1) {
+                  histLength1FakeNoCl[len - 4]->Fill(iLayer);
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+
+    auto canvas = new TCanvas("fc_canvas", "Fake clusters", 1600, 1000);
+    canvas->Divide(4, 2);
+    for (int iH{0}; iH < 4; ++iH) {
+      canvas->cd(iH + 1);
+      stackLength[iH]->Draw();
+      gPad->BuildLegend();
+    }
+    for (int iH{0}; iH < 4; ++iH) {
+      canvas->cd(iH + 5);
+      stackLength1Fake[iH]->Draw();
+      gPad->BuildLegend();
+    }
+    canvas->SaveAs("fakeClusters.png", "recreate");
+  }
 }
@@ -26,6 +26,7 @@ o2_add_library(ITStracking
                       src/ClusterLines.cxx
                       src/Vertexer.cxx
                       src/VertexerTraits.cxx
+                       src/Smoother.cxx
               PUBLIC_LINK_LIBRARIES O2::GPUCommon
                                     Microsoft.GSL::GSL
                                     O2::CommonConstants

@@ -0,0 +1,69 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+//
+/// \file Smoother.h
+/// \brief Class to handle Kalman smoothing for ITS tracking.
+///         Its instance stores the state of the track to the level we want to smooth to avoid multiple re-propagations when testing different clusters.
+///
+
+#include "ReconstructionDataFormats/Track.h"
+#include "DataFormatsITS/TrackITS.h"
+#include "DetectorsBase/Propagator.h"
+#include "ITStracking/ROframe.h"
+
+#if defined(CA_DEBUG) || defined(CA_STANDALONE_DEBUGGER)
+#include "ITStracking/StandaloneDebugger.h"
+#endif
+
+namespace o2
+{
+namespace its
+{
+
+template <unsigned int D>
+class Smoother
+{
+ public:
+  Smoother(TrackITSExt& track, int layer, const ROframe& event, float bZ, o2::base::PropagatorF::MatCorrType corr);
+  ~Smoother();
+
+  bool isValidInit() const
+  {
+    return mInitStatus;
+  }
+  bool testCluster(const int clusterId, const ROframe& event);
+  bool getSmoothedTrack();
+  float getChi2() const { return mBestChi2; }
+  float getLastChi2() const { return mLastChi2; }
+
+ private:
+  float computeSmoothedPredictedChi2(const o2::track::TrackParCov& outwTrack,
+                                     const o2::track::TrackParCov& inwTrack,
+                                     const std::array<float, 2>& cls,
+                                     const std::array<float, 3>& clCov);
+  bool smoothTrack();
+
+ private:
+  int mLayerToSmooth;                       // Layer to compute smoothing optimization
+  float mBz;                                // Magnetic field along Z
+  bool mInitStatus;                         // State after the initialization
+  o2::base::PropagatorF::MatCorrType mCorr; // Type of correction to use
+  TrackITSExt mInwardsTrack;                // outwards track: from innermost cluster to outermost
+  TrackITSExt mOutwardsTrack;               // inwards track: from outermost cluster to innermost
+  float mBestChi2;                          // Best value of local smoothed chi2
+  float mLastChi2 = 1e8;                    // Latest computed chi2
+
+#if defined(CA_DEBUG) || defined(CA_STANDALONE_DEBUGGER)
+  StandaloneDebugger* mDebugger;
+#endif
+};
+} // namespace its
+} // namespace o2
@@ -70,6 +70,10 @@ class Tracker
  float getBz() const;

  void clustersToTracks(std::function<void(std::string s)> = [](std::string s) { std::cout << s << std::endl; });
+  void setSmoothing(bool v) { mApplySmoothing = v; }
+  bool getSmoothing() const { return mApplySmoothing; }
+
+  std::vector<TrackITSExt>& getTracks();

  void setROFrame(std::uint32_t f) { mROFrame = f; }
  std::uint32_t getROFrame() const { return mROFrame; }
@@ -104,6 +108,7 @@ class Tracker
  std::vector<TrackingParameters> mTrkParams;

  bool mCUDA = false;
+  bool mApplySmoothing = false;
  o2::base::PropagatorImpl<float>::MatCorrType mCorrType = o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrLUT;
  float mBz = 5.f;
  std::uint32_t mROFrame = 0;