Experimantal STS-Pixel Digitizer and Hitfinder

core/data/sts/CbmStsHit.h

@@ -73,6 +73,14 @@ public:
      **/
   double GetDu() const { return fDu; }
 
+  /** @brief Error of coordinate across front-side strips
+     ** @value Coordinate error [cm]
+     **
+     ** Note that this error is defined only in the
+     ** local coordinate system of the sensor.
+     **/
+  void SetDu(Double_t du) { fDu = du; }
+
 
   /** @brief Error of coordinate across front-side strips
      ** @value Coordinate error [cm]
@@ -82,6 +90,14 @@ public:
      **/
   double GetDv() const { return fDv; }
 
+  /** @brief Error of coordinate across front-side strips
+     ** @value Coordinate error [cm]
+     **
+     ** Note that this error is defined only in the
+     ** local coordinate system of the sensor.
+     **/
+  void SetDv(Double_t dv) { fDv = dv; }
+
 
   /** Index of cluster at the front side
      ** @value  Front-side cluster index

reco/L1/CbmL1.h

@@ -232,7 +232,7 @@ public:
   L1Algo* algo {nullptr};  // for access to L1 Algorithm from L1::Instance
 
   TString fMuchDigiFile {};  // Much digitization file name
-  bool fUseHitErrors {false};
+  bool fUseHitErrors {true};
   bool fMissingHits {false};
   L1Algo::TrackingMode fTrackingMode {L1Algo::TrackingMode::kSts};
 

reco/L1/L1Algo/L1Algo.h

@@ -256,7 +256,7 @@ public:
   L1Vector<int> fStripToTrackB {"L1Algo::fStripToTrackB"};  // back strip to track pointers
 
   int fNThreads {0};
-  bool fUseHitErrors {0};
+  bool fUseHitErrors {true};
   bool fMissingHits {0};
   TrackingMode fTrackingMode {kSts};
 

reco/detectors/sts/CMakeLists.txt

@@ -17,6 +17,7 @@ CbmStsAlgoFindHitsOrtho.cxx
 CbmStsFindTracks.cxx
 CbmStsFindTracksEvents.cxx
 CbmStsRecoModule.cxx
+CbmRecoStsPixel.cxx
 CbmStsTrackFinderIdeal.cxx
 
 unpack/CbmStsUnpackAlgo.cxx

reco/detectors/sts/CbmRecoStsLinkDef.h

@@ -17,6 +17,9 @@
 #pragma link C++ class CbmStsFindTracksEvents + ;
 #pragma link C++ class CbmStsFindTracksQa + ;
 #pragma link C++ class CbmStsRecoModule + ;
+
+#pragma link C++ class CbmRecoStsPixel + ;
+
 #pragma link C++ class CbmStsTrackFinderIdeal + ;
 
 #pragma link C++ class CbmStsUnpackAlgo + ;

reco/detectors/sts/CbmRecoStsPixel.cxx

+/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer] */
+
+/** @file CbmRecoStsPixel.cxx
+ ** @author Sergey Gorbunov
+ ** @since 09.12.2021
+ **/
+
+#include "CbmRecoStsPixel.h"
+
+#include "CbmAddress.h"
+#include "CbmDigiManager.h"
+#include "CbmEvent.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCDataManager.h"
+#include "CbmStsDigi.h"
+#include "CbmStsModule.h"
+#include "CbmStsParSetModule.h"
+#include "CbmStsParSetSensor.h"
+#include "CbmStsParSetSensorCond.h"
+#include "CbmStsParSim.h"
+#include "CbmStsPoint.h"
+#include "CbmStsRecoModule.h"
+#include "CbmStsSetup.h"
+#include "CbmStsStation.h"
+
+#include <FairField.h>
+#include <FairRun.h>
+#include <FairRuntimeDb.h>
+
+#include <TClonesArray.h>
+#include <TGeoBBox.h>
+#include <TGeoPhysicalNode.h>
+#include <TRandom.h>
+
+#include <iomanip>
+
+using std::fixed;
+using std::left;
+using std::right;
+using std::setprecision;
+using std::setw;
+using std::stringstream;
+using std::vector;
+
+
+ClassImp(CbmRecoStsPixel)
+
+
+  // -----   Constructor   ---------------------------------------------------
+  CbmRecoStsPixel::CbmRecoStsPixel(ECbmRecoMode mode)
+  : FairTask("RecoStsPixel", 1)
+  , fMode(mode)
+{
+}
+// -------------------------------------------------------------------------
+
+
+// -----   Destructor   ----------------------------------------------------
+CbmRecoStsPixel::~CbmRecoStsPixel() {}
+// -------------------------------------------------------------------------
+
+
+// -----   End-of-run action   ---------------------------------------------
+void CbmRecoStsPixel::Finish() {}
+// -------------------------------------------------------------------------
+
+
+// -----   Initialisation   ------------------------------------------------
+InitStatus CbmRecoStsPixel::Init()
+{
+
+  // --- Something for the screen
+  std::cout << std::endl;
+  LOG(info) << "==========================================================";
+  LOG(info) << GetName() << ": Initialising ";
+
+  // --- Check IO-Manager
+  FairRootManager* ioman = FairRootManager::Instance();
+  assert(ioman);
+
+  // --- In event mode: get input array (CbmEvent)
+  if (fMode == kCbmRecoEvent) {
+    LOG(info) << GetName() << ": Using event-by-event mode";
+    fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
+    if (nullptr == fEvents) {
+      LOG(warn) << GetName() << ": Event mode selected but no event array found!";
+      return kFATAL;
+    }  //? Event branch not present
+  }    //? Event mode
+  else {
+    LOG(info) << GetName() << ": Using time-based mode";
+  }
+
+
+  // --- Digi Manager
+  fDigiManager = CbmDigiManager::Instance();
+  fDigiManager->Init();
+
+  // --- Check input array (StsDigis)
+  if (!fDigiManager->IsPresent(ECbmModuleId::kSts)) {
+    LOG(fatal) << GetName() << ": No StsDigi branch in input!";
+    return kERROR;
+  }
+
+  if (!fDigiManager->IsMatchPresent(ECbmModuleId::kSts)) {
+    LOG(error) << GetName() << " sts digi matches are not present";
+    return kERROR;
+  }
+
+  CbmMCDataManager* mcManager = (CbmMCDataManager*) ioman->GetObject("MCDataManager");
+  if (!mcManager) {
+    LOG(error) << GetName() << ": No CbmMCDataManager!";
+    return kERROR;
+  }
+
+  fStsPoints = mcManager->InitBranch("StsPoint");
+
+  if (!fStsPoints) {
+    LOG(fatal) << GetName() << ": No StsPoint branch in input!";
+    return kERROR;
+  }
+
+
+  // --- Register output array
+  fClusters = new TClonesArray("CbmStsCluster", 1);
+  ioman->Register("StsCluster", "Clusters in STS", fClusters, IsOutputBranchPersistent("StsCluster"));
+
+  // --- Register output array
+  fHits = new TClonesArray("CbmStsHit", 1);
+  ioman->Register("StsHit", "Hits in STS", fHits, IsOutputBranchPersistent("StsHit"));
+
+  // --- Simulation settings
+  assert(fParSim);
+  LOG(info) << GetName() << ": Sim settings " << fParSim->ToString();
+
+  // --- Module parameters
+  assert(fParSetModule);
+  LOG(info) << GetName() << ": Module parameters " << fParSetModule->ToString();
+
+  // --- Sensor parameters
+  assert(fParSetSensor);
+  LOG(info) << GetName() << ": Sensor parameters " << fParSetModule->ToString();
+
+  // --- Sensor conditions
+  assert(fParSetCond);
+  //assert(fParSetCond->IsSet());
+  LOG(info) << GetName() << ": Sensor conditions " << fParSetCond->ToString();
+
+  // --- Initialise STS setup
+  fSetup = CbmStsSetup::Instance();
+  if (!fSetup->IsInit()) { fSetup->Init(nullptr); }
+  if (!fSetup->IsModuleParsInit()) { fSetup->SetModuleParameters(fParSetModule); }
+  if (!fSetup->IsSensorParsInit()) { fSetup->SetSensorParameters(fParSetSensor); }
+  if (!fSetup->IsSensorCondInit()) { fSetup->SetSensorConditions(fParSetCond); }
+
+  // make sure that the STS digis were produced by the experimental Pixel digitiser
+  if (strcmp(fParSetSensor->getDescription(), "Experimental STS Pixels")) {
+    LOG(error) << GetName() << " STS digis must be produced by the CbmStsDigitizePixel digitizer";
+    return kERROR;
+  }
+
+  LOG(info) << GetName() << ": Initialisation successful.";
+  LOG(info) << "==========================================================";
+
+  return kSUCCESS;
+}
+// -------------------------------------------------------------------------
+
+
+// -----   Task execution   ------------------------------------------------
+void CbmRecoStsPixel::Exec(Option_t*)
+{
+
+  // --- Clear hit output array
+  fHits->Delete();
+
+  // --- Reset cluster output array
+  fClusters->Delete();
+
+  if (fMode == kCbmRecoTimeslice) {
+    // --- Time-slice mode: process entire array
+    ProcessData(nullptr);
+  }
+  else {
+    // --- Event mode: loop over events
+    assert(fEvents);
+    for (Int_t iEvent = 0; iEvent < fEvents->GetEntriesFast(); iEvent++) {
+      CbmEvent* event = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));
+      assert(event);
+      ProcessData(event);
+    }  //# events
+  }    //? event mode
+}
+// -------------------------------------------------------------------------
+
+
+// -----   Process one time slice or event   -------------------------------
+void CbmRecoStsPixel::ProcessData(CbmEvent* event)
+{
+  if (!fDigiManager->IsMatchPresent(ECbmModuleId::kSts)) {
+    LOG(error) << GetName() << ": No StsDigi branch in input!";
+    return;
+  }
+
+  Int_t nDigis = (event ? event->GetNofData(ECbmDataType::kStsDigi) : fDigiManager->GetNofDigis(ECbmModuleId::kSts));
+  int nHits    = 0;
+
+  for (Int_t iDigi = 0; iDigi < nDigis; iDigi++) {
+
+    UInt_t digiIndex       = (event ? event->GetIndex(ECbmDataType::kStsDigi, iDigi) : iDigi);
+    const CbmStsDigi* digi = fDigiManager->Get<const CbmStsDigi>(digiIndex);
+    assert(digi);
+
+    // Check system ID. There are pulser digis in which will be ignored here.
+    Int_t systemId = CbmAddress::GetSystemId(digi->GetAddress());
+    if (systemId != ToIntegralType(ECbmModuleId::kSts)) { continue; }
+
+    const CbmMatch* match = fDigiManager->GetMatch(ECbmModuleId::kSts, digiIndex);
+    assert(match);
+
+    // make sure that the digi was produced by CbmStsDigitizePixel task
+    if (digi->GetChannel() != 0 || match->GetNofLinks() != 1) {
+      LOG(error) << GetName() << " sts digis were not produced by CbmStsDigitizePixel task";
+      break;
+    }
+
+    // the digi was produced by one MC point
+
+    const CbmLink& link = match->GetLink(0);
+    CbmStsPoint* pt     = dynamic_cast<CbmStsPoint*>(fStsPoints->Get(link.GetFile(), link.GetEntry(), link.GetIndex()));
+
+    // create one cluster and one hit per digi
+
+    UInt_t iCluster        = fClusters->GetEntriesFast();
+    CbmStsCluster* cluster = new ((*fClusters)[iCluster]) CbmStsCluster;
+    cluster->AddDigi(iDigi);
+    if (event) event->AddData(ECbmDataType::kStsCluster, iCluster);
+
+    UInt_t iHit    = fHits->GetEntriesFast();
+    CbmStsHit* hit = new ((*fHits)[iHit]) CbmStsHit;
+
+    hit->SetAddress(digi->GetAddress());
+    hit->SetFrontClusterId(iCluster);
+    hit->SetBackClusterId(iCluster);
+
+    int ista = CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress());
+    if (ista < 0 || ista >= CbmStsSetup::Instance()->GetNofStations()) {
+      LOG(error) << "wrong Sts station number " << ista;
+      break;
+    }
+
+    CbmStsStation* station = CbmStsSetup::Instance()->GetStation(ista);
+    double staZ            = station->GetZ();
+
+    if (fabs(pt->GetZ() - staZ) > 1.) {
+      LOG(error) << "Sts point Z " << pt->GetZ() << " is far from the station Z " << staZ;
+      break;
+    }
+
+    double dx = pt->GetXOut() - pt->GetXIn();
+    double dy = pt->GetYOut() - pt->GetYIn();
+    double dz = pt->GetZOut() - pt->GetZIn();
+    if (fabs(dz) > 1.e-2) {
+      hit->SetX(pt->GetXIn() + dx * (staZ - pt->GetZIn()) / dz);
+      hit->SetY(pt->GetYIn() + dy * (staZ - pt->GetZIn()) / dz);
+      hit->SetZ(staZ);
+    }
+    else {
+      hit->SetX(pt->GetXIn());
+      hit->SetY(pt->GetYIn());
+      hit->SetZ(pt->GetZIn());
+    }
+
+    Double_t pitchX = station->GetSensorPitch(0);
+    Double_t pitchY = station->GetSensorPitch(1);
+
+    assert(pitchX > 1.e-5);
+    assert(pitchY > 1.e-5);
+
+    hit->SetX((floor(hit->GetX() / pitchX) + 0.5) * pitchX);
+    hit->SetY((floor(hit->GetY() / pitchY) + 0.5) * pitchY);
+
+    hit->SetDx(pitchX / sqrt(12.));
+    hit->SetDy(pitchY / sqrt(12.));
+    hit->SetDxy(0.);
+    hit->SetDu(hit->GetDx());
+    hit->SetDv(hit->GetDy());
+
+    hit->SetTime(digi->GetTime());
+    const CbmStsParModule& module = fParSetModule->GetParModule(digi->GetAddress());
+    const CbmStsParAsic& asic     = module.GetParAsic(0);
+    hit->SetTimeError(asic.GetTimeResol());
+
+    if (event) event->AddData(ECbmDataType::kStsHit, iHit);
+    nHits++;
+  }  // digis
+
+  LOG(info) << GetName() << " Processed " << nDigis << " digis, created " << nHits << " hits";
+}
+// -------------------------------------------------------------------------
+
+
+// -----   Connect parameter container   -----------------------------------
+void CbmRecoStsPixel::SetParContainers()
+{
+  FairRuntimeDb* db = FairRun::Instance()->GetRuntimeDb();
+  fParSim           = dynamic_cast<CbmStsParSim*>(db->getContainer("CbmStsParSim"));
+  fParSetModule     = dynamic_cast<CbmStsParSetModule*>(db->getContainer("CbmStsParSetModule"));
+  fParSetSensor     = dynamic_cast<CbmStsParSetSensor*>(db->getContainer("CbmStsParSetSensor"));
+  fParSetCond       = dynamic_cast<CbmStsParSetSensorCond*>(db->getContainer("CbmStsParSetSensorCond"));
+}
+// -------------------------------------------------------------------------

reco/detectors/sts/CbmRecoStsPixel.h

+/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer] */
+
+/** @file CbmRecoStsPixel.h
+ ** @author Sergey Gorbunov
+ ** @since 09.12.2021
+ **/
+
+
+#ifndef CbmRecoStsPixel_H
+#define CbmRecoStsPixel_H 1
+
+#include "CbmRecoSts.h"
+
+#include <FairTask.h>
+
+#include <TStopwatch.h>
+
+
+class CbmDigiManager;
+class CbmEvent;
+class CbmStsElement;
+class CbmStsParAsic;
+class CbmStsParModule;
+class CbmStsParSensor;
+class CbmStsParSensorCond;
+class CbmStsParSetModule;
+class CbmStsParSetSensor;
+class CbmStsParSetSensorCond;
+class CbmStsParSim;
+class CbmStsRecoModule;
+class CbmStsSensor;
+class CbmStsSetup;
+class CbmMCDataArray;
+
+/** @class CbmRecoStsPixel
+ ** @brief Task class for local reconstruction in the STS Pixel detector
+ ** @author Sergey Gorbunov
+ ** @since 09.12.2021
+ ** @date 09.12.2021
+ **
+ ** Local reconstruction for the experimental STS Pixel detector
+ **
+ ** The STS digis must be produced by the CbmStsDigitizePixel task.
+ ** To run it, replace CbmRecoSts task by CbmRecoStsPixel task in the run_reco.C macro
+ **
+ **/
+class CbmRecoStsPixel : public FairTask {
+
+public:
+  /** @brief Constructor **/
+  CbmRecoStsPixel(ECbmRecoMode mode = kCbmRecoTimeslice);
+
+
+  /** @brief Copy constructor (disabled) **/
+  CbmRecoStsPixel(const CbmRecoStsPixel&) = delete;
+
+
+  /** @brief Assignment operator (disabled) **/
+  CbmRecoStsPixel operator=(const CbmRecoStsPixel&) = delete;
+
+
+  /** @brief Destructor  **/
+  ~CbmRecoStsPixel();
+
+
+  /** @brief Task execution **/
+  void Exec(Option_t* opt);
+
+
+  /** @brief End-of-run action **/
+  void Finish();
+
+
+  /** @brief Initialisation **/
+  InitStatus Init();
+
+
+  /** @brief Set event-by-event mode
+     ** @param choice  If true, event-by-event mode is used
+     **
+     ** In the event-by-event mode, the event objects in the input tree
+     ** are used, and events are processed one after the other.
+     ** An event builder has to be run before, creating the event objects.
+     ** By default, time-slice mode is applied.
+     **
+     ** Alternative to using SetMode.
+     **/
+  void SetEventMode(Bool_t choice = kTRUE) { fMode = (choice ? kCbmRecoEvent : kCbmRecoTimeslice); }
+
+
+  /** @brief Set execution mode
+     ** @param mode  Time-slice or event
+     **
+     ** In the time-slice mode, the entire time-slice (input arrays)
+     ** will be processed. In the event mode, events read from the event
+     ** branch are processed one after the other.
+     **/
+  void SetMode(ECbmRecoMode mode) { fMode = mode; }
+
+
+  /** @brief Define the needed parameter containers **/
+  void SetParContainers();
+
+
+private:
+  /** @brief Process one time slice or event
+     ** @param event  Pointer to CbmEvent object
+     **
+     ** If a null event pointer is given, the entire input array is processed.
+     **/
+  void ProcessData(CbmEvent* event = nullptr);
+
+private:
+  // --- I/O
+  TClonesArray* fEvents        = nullptr;  //! Input array of events
+  CbmDigiManager* fDigiManager = nullptr;  //! Interface to digi branch
+  CbmMCDataArray* fStsPoints {nullptr};
+  TClonesArray* fClusters = nullptr;  //! Output cluster array
+  TClonesArray* fHits     = nullptr;  //! Output hit array
+
+  // --- Setup and parameters
+  CbmStsSetup* fSetup                 = nullptr;  //! Instance of STS setup
+  CbmStsParSim* fParSim               = nullptr;  ///< Simulation settings
+  CbmStsParSetModule* fParSetModule   = nullptr;  ///< Module parameters
+  CbmStsParSetSensor* fParSetSensor   = nullptr;  ///< Sensor parameters
+  CbmStsParSetSensorCond* fParSetCond = nullptr;  ///< Sensor conditions
+
+
+  // --- Settings
+  ECbmRecoMode fMode = kCbmRecoEvent;  ///< Time-slice or event
+
+  ClassDef(CbmRecoStsPixel, 0);
+};
+
+#endif

sim/detectors/sts/CMakeLists.txt

@@ -20,6 +20,8 @@ CbmStsSimSensorDssdOrtho.cxx
 CbmStsSimSensorDssdStereo.cxx
 CbmStsSimSensorFactory.cxx
 
+CbmStsDigitizePixel.cxx
+
 qa/CbmStsDigitizeQa.cxx
 qa/CbmStsDigitizeQaReport.cxx
 )
@@ -103,6 +105,7 @@ FILES CbmStsDigitize.h
 	  CbmStsSimSensorDssdOrtho.h
 	  CbmStsSimSensorFactory.h
 	  CbmStsTrackStatus.h
+	  CbmStsDigitizePixel.h
 DESTINATION include
 )
 # ---------------------------------------------------------

sim/detectors/sts/CbmStsDigitizePixel.cxx

+/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer] */
+
+/** @file CbmStsDigitizePixel.cxx
+ ** @author Sergey Gorbunov
+ ** @date 09.12.2021
+ **/
+
+// Include class header
+#include "CbmStsDigitizePixel.h"
+
+// Includes from C++
+#include <cassert>
+#include <cstring>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+
+// Includes from ROOT
+#include "TClonesArray.h"
+#include "TGeoBBox.h"
+#include "TGeoMatrix.h"
+#include "TGeoPhysicalNode.h"
+#include "TGeoVolume.h"
+#include <TMCProcess.h>
+
+// Includes from FairRoot
+#include "FairEventHeader.h"
+#include "FairField.h"
+#include "FairLink.h"
+#include "FairMCEventHeader.h"
+#include "FairMCPoint.h"
+#include "FairRunAna.h"
+#include "FairRunSim.h"
+#include "FairRuntimeDb.h"
+#include <Logger.h>
+
+// Includes from CbmRoot
+#include "CbmStsDigi.h"
+#include "CbmStsParSetModule.h"
+#include "CbmStsPoint.h"
+#include "CbmStsSetup.h"
+
+// Includes from STS
+#include "CbmStsModule.h"
+#include "CbmStsParAsic.h"
+#include "CbmStsParModule.h"
+#include "CbmStsParSensor.h"
+#include "CbmStsParSensorCond.h"
+#include "CbmStsParSetModule.h"
+#include "CbmStsParSetSensor.h"
+#include "CbmStsParSetSensorCond.h"
+#include "CbmStsParSim.h"
+#include "CbmStsPhysics.h"
+#include "CbmStsSensor.h"
+#include "CbmStsSetup.h"
+#include "CbmStsSimSensorFactory.h"
+
+
+using std::fixed;
+using std::left;
+using std::right;
+using std::setprecision;
+using std::setw;
+using std::string;
+using std::stringstream;
+
+using namespace CbmSts;
+
+ClassImp(CbmStsDigitizePixel);