diff --git a/macro/L1/run_reco_L1global.C b/macro/L1/run_reco_L1global.C
new file mode 100644
index 0000000000000000000000000000000000000000..7346ec510d0cf5c81786936cce8edfa8806ec788
--- /dev/null
+++ b/macro/L1/run_reco_L1global.C
@@ -0,0 +1,520 @@
+/* Copyright (C) 2020-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Volker Friese [committer], Dominik Smith */
+
+/** @file run_reco_L1global.C
+ ** @author Sergey Gorbunov
+ ** @since 1 Juni 2022
+ **/
+
+
+// --- Includes needed for IDE
+#include <RtypesCore.h>
+#if !defined(__CLING__)
+#include "CbmBuildEventsFromTracksReal.h"
+#include "CbmBuildEventsIdeal.h"
+#include "CbmBuildEventsQa.h"
+#include "CbmDefs.h"
+#include "CbmFindPrimaryVertex.h"
+#include "CbmKF.h"
+#include "CbmL1.h"
+#include "CbmL1StsTrackFinder.h"
+#include "CbmLitFindGlobalTracks.h"
+#include "CbmMCDataManager.h"
+#include "CbmMatchRecoToMC.h"
+#include "CbmMuchFindHitsGem.h"
+#include "CbmMvdClusterfinder.h"
+#include "CbmMvdHitfinder.h"
+#include "CbmPVFinderKF.h"
+#include "CbmPrimaryVertexFinder.h"
+#include "CbmPsdHitProducer.h"
+#include "CbmRecoSts.h"
+#include "CbmRichHitProducer.h"
+#include "CbmRichReconstruction.h"
+#include "CbmSetup.h"
+#include "CbmStsFindTracks.h"
+#include "CbmStsFindTracksEvents.h"
+#include "CbmStsTrackFinder.h"
+#include "CbmTaskBuildRawEvents.h"
+#include "CbmTofSimpClusterizer.h"
+#include "CbmTrdClusterFinder.h"
+#include "CbmTrdHitProducer.h"
+
+#include <FairFileSource.h>
+#include <FairMonitor.h>
+#include <FairParAsciiFileIo.h>
+#include <FairParRootFileIo.h>
+#include <FairRunAna.h>
+#include <FairRuntimeDb.h>
+#include <FairSystemInfo.h>
+
+#include <TStopwatch.h>
+#endif
+
+
+/** @brief Macro for CBM reconstruction
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @since 14 November 2020
+ ** @param input Name of input file (w/o extension .raw.root)
+ ** @param nTimeSlices Number of time-slices to process
+ ** @param firstTimeSlice First time-slice (entry) to be processed
+ ** @param output Name of output file (w/o extension .rec.root)
+ ** @param sEvBuildRaw Option for raw event building
+ ** @param setup Name of predefined geometry setup
+ ** @param paramFile Parameter ROOT file (w/o extension .par.root)
+ ** @param debugWithMC Option to provide the trackfinder with MC information
+ **
+ ** This macro is a copy of macro/reco/run_reco.C
+ ** with a test version of L1 global tracker
+ **
+ **/
+void run_reco_L1global(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice = 0, TString output = "",
+ TString sEvBuildRaw = "", TString setup = "sis100_electron", TString paramFile = "",
+ Bool_t debugWithMC = false)
+{
+
+ // how to run it:
+ // root -l -q $VMCWORKDIR/macro/L1/run_reco_L1global.C'("trd2d", -1, 0, "", "", "trd2d", "", true)'
+
+ // ========================================================================
+ // Adjust this part according to your requirements
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "INFO";
+ TString logVerbosity = "LOW";
+ // ------------------------------------------------------------------------
+
+ // ----- Environment --------------------------------------------------
+ TString myName = "run_reco"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+
+ // ----- In- and output file names ------------------------------------
+ if (input.IsNull()) input = "test";
+ TString rawFile = input + ".raw.root";
+ TString traFile = input + ".tra.root";
+ if (output.IsNull()) output = input;
+ TString outFile = output + ".reco.root";
+ TString monFile = output + ".moni_reco.root";
+ if (paramFile.IsNull()) paramFile = input;
+ TString parFile = paramFile + ".par.root";
+ std::cout << "Inputfile " << rawFile << std::endl;
+ std::cout << "Outfile " << outFile << std::endl;
+ std::cout << "Parfile " << parFile << std::endl;
+
+ // ----- Load the geometry setup -------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl;
+ CbmSetup* geo = CbmSetup::Instance();
+ geo->LoadSetup(setup);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Some global switches -----------------------------------------
+ Bool_t eventBased = !sEvBuildRaw.IsNull();
+ Bool_t useMvd = geo->IsActive(ECbmModuleId::kMvd);
+ Bool_t useSts = geo->IsActive(ECbmModuleId::kSts);
+ Bool_t useRich = geo->IsActive(ECbmModuleId::kRich);
+ Bool_t useMuch = geo->IsActive(ECbmModuleId::kMuch);
+ Bool_t useTrd = geo->IsActive(ECbmModuleId::kTrd);
+ Bool_t useTof = geo->IsActive(ECbmModuleId::kTof);
+ Bool_t usePsd = geo->IsActive(ECbmModuleId::kPsd);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Parameter files as input to the runtime database -------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Defining parameter files " << std::endl;
+ TList* parFileList = new TList();
+ TString geoTag;
+
+ // - TRD digitisation parameters
+ if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTrd, geoTag)) {
+ const Char_t* npar[4] = {"asic", "digi", "gas", "gain"};
+ TObjString* trdParFile(NULL);
+ for (Int_t i(0); i < 4; i++) {
+ trdParFile = new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + "." + npar[i] + ".par");
+ parFileList->Add(trdParFile);
+ std::cout << "-I- " << myName << ": Using parameter file " << trdParFile->GetString() << std::endl;
+ }
+ }
+
+ // - TOF digitisation parameters
+ if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTof, geoTag)) {
+ TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
+ parFileList->Add(tofBdfFile);
+ std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+ // In general, the following parts need not be touched
+ // ========================================================================
+
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+
+ // ----- FairRunAna ---------------------------------------------------
+ FairRunAna* run = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(rawFile);
+ if (debugWithMC) { inputSource->AddFriend(traFile); }
+ run->SetSource(inputSource);
+ run->SetOutputFile(outFile);
+ run->SetGenerateRunInfo(kTRUE);
+ FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monFile);
+ // ------------------------------------------------------------------------
+
+ // ----- MCDataManager -----------------------------------
+ if (debugWithMC) {
+ CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 0);
+ mcManager->AddFile(traFile);
+ run->AddTask(mcManager);
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- Logger settings ----------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Raw event building from digis --------------------------------
+ if (eventBased) {
+ if (sEvBuildRaw.EqualTo("Ideal", TString::ECaseCompare::kIgnoreCase)) {
+ FairTask* evBuildRaw = new CbmBuildEventsIdeal();
+ run->AddTask(evBuildRaw);
+ std::cout << "-I- " << myName << ": Added task " << evBuildRaw->GetName() << std::endl;
+ } //? Ideal raw event building
+ else if (sEvBuildRaw.EqualTo("Real", TString::ECaseCompare::kIgnoreCase)) {
+ CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents();
+
+ //Choose between NoOverlap, MergeOverlap, AllowOverlap
+ evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap);
+
+ // Remove detectors where digis not found
+ if (!useRich) evBuildRaw->RemoveDetector(kRawEventBuilderDetRich);
+ if (!useMuch) evBuildRaw->RemoveDetector(kRawEventBuilderDetMuch);
+ if (!usePsd) evBuildRaw->RemoveDetector(kRawEventBuilderDetPsd);
+ if (!useTof) evBuildRaw->RemoveDetector(kRawEventBuilderDetTof);
+ if (!useTrd) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd);
+ if (!useSts) {
+ std::cerr << "-E- " << myName << ": Sts must be present for raw event "
+ << "building using ``Real2019'' option. Terminating macro." << std::endl;
+ return;
+ }
+ // Set STS as reference detector
+ evBuildRaw->SetReferenceDetector(kRawEventBuilderDetSts);
+
+ // Use sliding window seed builder with STS
+ //evBuildRaw->SetReferenceDetector(kRawEventBuilderDetUndef);
+ //evBuildRaw->AddSeedTimeFillerToList(kRawEventBuilderDetSts);
+ //evBuildRaw->SetSlidingWindowSeedFinder(1000, 500, 500);
+ //evBuildRaw->SetSeedFinderQa(true); // optional QA information for seed finder
+
+ evBuildRaw->SetTsParameters(0.0, 1.e7, 0.0);
+
+ // Use CbmMuchDigi instead of CbmMuchBeamtimeDigi
+ evBuildRaw->ChangeMuchBeamtimeDigiFlag(kFALSE);
+
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, 1000);
+ evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1);
+ evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, -500, 500);
+
+ run->AddTask(evBuildRaw);
+ std::cout << "-I- " << myName << ": Added task " << evBuildRaw->GetName() << std::endl;
+ } //? Real raw event building
+ else {
+ std::cerr << "-E- " << myName << ": Unknown option " << sEvBuildRaw
+ << " for raw event building! Terminating macro execution." << std::endl;
+ return;
+ }
+ } //? event-based reco
+ // ------------------------------------------------------------------------
+
+
+ // ----------- QA for raw event builder -----------------------------------
+ if (eventBased && debugWithMC) {
+ CbmBuildEventsQa* evBuildQA = new CbmBuildEventsQa();
+ run->AddTask(evBuildQA);
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in MVD ----------------------------------
+ if (useMvd) {
+
+ CbmMvdClusterfinder* mvdCluster = new CbmMvdClusterfinder("MVD Cluster Finder", 0, 0);
+ run->AddTask(mvdCluster);
+ std::cout << "-I- : Added task " << mvdCluster->GetName() << std::endl;
+
+ CbmMvdHitfinder* mvdHit = new CbmMvdHitfinder("MVD Hit Finder", 0, 0);
+ mvdHit->UseClusterfinder(kTRUE);
+ run->AddTask(mvdHit);
+ std::cout << "-I- " << myName << ": Added task " << mvdHit->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in STS ----------------------------------
+ if (useSts) {
+ CbmRecoSts* stsReco = new CbmRecoSts(ECbmRecoMode::kCbmRecoTimeslice);
+ if (eventBased) stsReco->SetMode(ECbmRecoMode::kCbmRecoEvent);
+ run->AddTask(stsReco);
+ std::cout << "-I- " << myName << ": Added task " << stsReco->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in RICH ---------------------------------
+ if (useRich) {
+ CbmRichHitProducer* richHitProd = new CbmRichHitProducer();
+ run->AddTask(richHitProd);
+ std::cout << "-I- " << myName << ": Added task " << richHitProd->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in MUCH ---------------------------------
+ if (useMuch) {
+
+ // --- Parameter file name
+ TString geoTag;
+ geo->GetGeoTag(ECbmModuleId::kMuch, geoTag);
+ Int_t muchFlag = (geoTag.Contains("mcbm") ? 1 : 0);
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile += "/parameters/much/much_" + geoTag(0, 4) + "_digi_sector.root";
+ std::cout << "Using parameter file " << parFile << std::endl;
+
+ // --- Hit finder for GEMs
+ FairTask* muchReco = new CbmMuchFindHitsGem(parFile.Data(), muchFlag);
+ run->AddTask(muchReco);
+ std::cout << "-I- " << myName << ": Added task " << muchReco->GetName() << " with parameter file " << parFile
+ << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in TRD ----------------------------------
+ if (useTrd) {
+
+ Double_t triggerThreshold = 0.5e-6; // SIS100
+ CbmTrdClusterFinder* trdCluster = new CbmTrdClusterFinder();
+ if (eventBased) trdCluster->SetTimeBased(kFALSE);
+ else
+ trdCluster->SetTimeBased(kTRUE);
+ trdCluster->SetNeighbourEnable(true, false);
+ trdCluster->SetMinimumChargeTH(triggerThreshold);
+ trdCluster->SetRowMerger(true);
+
+ // Uncomment if you want to use all available digis.
+ // In that case clusters hits will not be added to the CbmEvent
+ // trdCluster->SetUseOnlyEventDigis(kFALSE);
+
+ run->AddTask(trdCluster);
+ std::cout << "-I- " << myName << ": Added task " << trdCluster->GetName() << std::endl;
+
+ CbmTrdHitProducer* trdHit = new CbmTrdHitProducer();
+ run->AddTask(trdHit);
+ std::cout << "-I- " << myName << ": Added task " << trdHit->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in TOF ----------------------------------
+ if (useTof) {
+ CbmTofSimpClusterizer* tofCluster = new CbmTofSimpClusterizer("TOF Simple Clusterizer", 0);
+ tofCluster->SetOutputBranchPersistent("TofHit", kTRUE);
+ tofCluster->SetOutputBranchPersistent("TofDigiMatch", kTRUE);
+ run->AddTask(tofCluster);
+ std::cout << "-I- " << myName << ": Added task " << tofCluster->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in PSD ----------------------------------
+ if (usePsd) {
+ CbmPsdHitProducer* psdHit = new CbmPsdHitProducer();
+ run->AddTask(psdHit);
+ std::cout << "-I- " << myName << ": Added task " << psdHit->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+ if (debugWithMC) {
+ CbmMatchRecoToMC* match1 = new CbmMatchRecoToMC();
+ run->AddTask(match1);
+ }
+
+ // ----- Track finding in STS (+ MVD) --------------------------------
+ if (useMvd || useSts) {
+ CbmKF* kalman = new CbmKF();
+ run->AddTask(kalman);
+ CbmL1* l1 = 0;
+ if (debugWithMC) { l1 = new CbmL1("L1", 2, 2); }
+ else {
+ l1 = new CbmL1("L1", 0);
+ }
+ l1->SetGlobalMode();
+ l1->SetUseHitErrors(true);
+ l1->SetUseMcHit(0, 0, 1, 0);
+
+ // --- Material budget file names
+ TString mvdGeoTag;
+ if (geo->GetGeoTag(ECbmModuleId::kMvd, mvdGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile += "/parameters/mvd/mvd_matbudget_" + mvdGeoTag + ".root";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetMvdMaterialBudgetFileName(parFile.Data());
+ }
+ TString stsGeoTag;
+ if (geo->GetGeoTag(ECbmModuleId::kSts, stsGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile += "/parameters/sts/sts_matbudget_" + stsGeoTag + ".root";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetStsMaterialBudgetFileName(parFile.Data());
+ }
+
+ TString trdGeoTag;
+ if (geo->GetGeoTag(ECbmModuleId::kTrd, trdGeoTag)) {
+ TString parFile = gSystem->Getenv("VMCWORKDIR");
+ parFile = parFile + "/parameters/trd/trd_matbudget_" + trdGeoTag + ".root ";
+ std::cout << "Using material budget file " << parFile << std::endl;
+ l1->SetTrdMaterialBudgetFileName(parFile.Data());
+ }
+
+ run->AddTask(l1);
+ std::cout << "-I- " << myName << ": Added task " << l1->GetName() << std::endl;
+
+ CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder();
+ FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder);
+ run->AddTask(globalFindTracks);
+ std::cout << "-I- " << myName << ": Added task " << globalFindTracks->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ==== From here on, the time-based and the event-based reconstruction
+ // ==== chains differ, since time-based version of primary vertex finding
+ // ==== and global tracking are not yet available. For time-based
+ // ==== reconstruction, a track-based event finder is used; no global
+ // ==== tracks are produced.
+
+ if (eventBased) {
+
+ // ----- Primary vertex finding -------------------------------------
+ CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF();
+ CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
+ run->AddTask(findVertex);
+ std::cout << "-I- " << myName << ": Added task " << findVertex->GetName() << std::endl;
+ // ----------------------------------------------------------------------
+
+
+ // --- Global track finding -----------------------------------------
+ CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks();
+ finder->SetTrackingType("branch");
+ finder->SetMergerType("nearest_hit");
+ run->AddTask(finder);
+ std::cout << "-I- : Added task " << finder->GetName() << std::endl;
+ // ----------------------------------------------------------------------
+
+ // --- Particle Id in TRD -----------------------------------------
+ if (useTrd) {
+ CbmTrdSetTracksPidLike* trdLI = new CbmTrdSetTracksPidLike("TRDLikelihood", "TRDLikelihood");
+ trdLI->SetUseMCInfo(kTRUE);
+ trdLI->SetUseMomDependence(kTRUE);
+ run->AddTask(trdLI);
+ std::cout << "-I- : Added task " << trdLI->GetName() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- RICH reconstruction ----------------------------------------
+ if (useRich) {
+ CbmRichReconstruction* richReco = new CbmRichReconstruction();
+ run->AddTask(richReco);
+ std::cout << "-I- : Added task " << richReco->GetName() << std::endl;
+ }
+ // ----------------------------------------------------------------------
+
+ } //? event-based reco
+
+ else {
+
+ // -----Â Event building from STS tracks -----------------------------
+ run->AddTask(new CbmBuildEventsFromTracksReal());
+ // ----------------------------------------------------------------------
+
+ } //? time-based reco
+
+
+ // ----- Parameter database --------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Set runtime DB" << std::endl;
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ FairParRootFileIo* parIo1 = new FairParRootFileIo();
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ parIo1->open(parFile.Data(), "UPDATE");
+ rtdb->setFirstInput(parIo1);
+ if (!parFileList->IsEmpty()) {
+ parIo2->open(parFileList, "in");
+ rtdb->setSecondInput(parIo2);
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Run initialisation -------------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+ rtdb->setOutput(parIo1);
+ rtdb->saveOutput();
+ rtdb->print();
+ // ------------------------------------------------------------------------
+
+
+ // ----- Register light ions (d, t, He3, He4) -------------------------
+ std::cout << std::endl;
+ TString registerLightIonsMacro = gSystem->Getenv("VMCWORKDIR");
+ registerLightIonsMacro += "/macro/KF/registerLightIons.C";
+ std::cout << "Loading macro " << registerLightIonsMacro << std::endl;
+ gROOT->LoadMacro(registerLightIonsMacro);
+ gROOT->ProcessLine("registerLightIons()");
+ // ------------------------------------------------------------------------
+
+ // ----- Start run ----------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(firstTimeSlice, nTimeSlices);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Finish -------------------------------------------------------
+ timer.Stop();
+ FairMonitor::GetMonitor()->Print();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished successfully." << std::endl;
+ std::cout << "Output file is " << outFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << std::endl;
+ FairSystemInfo sysInfo;
+ Float_t maxMemory = sysInfo.GetMaxMemory();
+ std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ std::cout << maxMemory;
+ std::cout << "</DartMeasurement>" << std::endl;
+ Float_t cpuUsage = ctime / rtime;
+ std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ std::cout << cpuUsage;
+ std::cout << "</DartMeasurement>" << std::endl;
+ // ------------------------------------------------------------------------
+
+
+ // ----- This is to prevent a malloc error when exiting ROOT ----------
+ // The source of the error is unknown. Related to TGeoManager.
+ RemoveGeoManager();
+ // ------------------------------------------------------------------------
+
+} // End of main macro function
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 412a71af9b06d87d2519e0db31ed65dc8a214c8d..e7cea72a07e16c174fce579874396b2a7bd6c638 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -118,30 +118,11 @@ CbmL1::~CbmL1()
void CbmL1::CheckDetectorPresence()
{
- Bool_t IsMuch = CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch);
- Bool_t IsTrd = CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd);
- Bool_t IsTof = CbmSetup::Instance()->IsActive(ECbmModuleId::kTof);
- //Bool_t IsSts = CbmSetup::Instance()->IsActive(ECbmModuleId::kSts);
- Bool_t IsMvd = CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd);
-
- /*
- TObjArray* topNodes = gGeoManager->GetTopNode()->GetNodes();
-
- for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast(); iTopNode++) {
- TGeoNode* topNode = static_cast<TGeoNode*>(topNodes->At(iTopNode));
-
- if (TString(topNode->GetName()).Contains("much")) IsMuch = 1;
- if (TString(topNode->GetName()).Contains("trd")) IsTrd = 1;
- if (TString(topNode->GetName()).Contains("tof")) IsTof = 1;
- //if (TString(topNode->GetName()).Contains("sts")) IsSts = 1;
- if (TString(topNode->GetName()).Contains("mvd")) IsMvd = 1;
- }
- */
-
- fUseMUCH = (fUseMUCH && IsMuch);
- fUseTRD = fUseTRD && IsTrd;
- fUseTOF = fUseTOF && IsTof;
- fUseMVD = fUseMVD && IsMvd;
+ fUseMVD = fUseMVD && CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd);
+ fUseSTS = fUseSTS && CbmSetup::Instance()->IsActive(ECbmModuleId::kSts);
+ fUseMUCH = fUseMUCH && CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch);
+ fUseTRD = fUseTRD && CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd);
+ fUseTOF = fUseTOF && CbmSetup::Instance()->IsActive(ECbmModuleId::kTof);
}
@@ -207,35 +188,35 @@ InitStatus CbmL1::Init()
#endif
}
- FairRootManager* fManger = FairRootManager::Instance();
-
- FairRunAna* Run = FairRunAna::Instance();
- {
- fUseMVD = 1;
- CbmStsFindTracks* FindTask = L1_DYNAMIC_CAST<CbmStsFindTracks*>(Run->GetTask("STSFindTracks"));
- if (FindTask) fUseMVD = FindTask->MvdUsage();
- // TODO: include of CbmSetup.h creates problems on Mac
- // if (!CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd)) { fUseMVD = false; }
- // N Mvd stations is read from the KF material
- if (CbmKF::Instance()->vMvdMaterial.size() == 0) { fUseMVD = false; }
- }
-
fHistoDir = gROOT->mkdir("L1");
-
// turn on reconstruction in sub-detectors
+ fUseMVD = 1;
+ fUseSTS = 1;
fUseMUCH = 0;
fUseTRD = 0;
fUseTOF = 0;
- if (fTrackingMode == L1Algo::TrackingMode::kMcbm) {
+ FairRootManager* fairManager = FairRootManager::Instance();
+ {
+ CbmStsFindTracks* findTask = L1_DYNAMIC_CAST<CbmStsFindTracks*>(FairRunAna::Instance()->GetTask("STSFindTracks"));
+ if (findTask) fUseMVD = findTask->MvdUsage();
+ if (CbmKF::Instance()->vMvdMaterial.size() == 0) { fUseMVD = false; }
+ }
+
+
+ if (L1Algo::TrackingMode::kMcbm == fTrackingMode) {
fUseMUCH = 1;
fUseTRD = 1;
fUseTOF = 1;
}
- if (fTrackingMode == L1Algo::TrackingMode::kGlobal) {
+
+ if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
+ //SGtrd2D!!
+ fUseMVD = 0;
+ fUseSTS = 0;
fUseMUCH = 0;
fUseTRD = 1;
fUseTOF = 0;
@@ -262,7 +243,7 @@ InitStatus CbmL1::Init()
if (!fLegacyEventMode) { // Time-slice mode selected
LOG(info) << GetName() << ": running in time-slice mode.";
fTimeSlice = NULL;
- fTimeSlice = (CbmTimeSlice*) fManger->GetObject("TimeSlice.");
+ fTimeSlice = (CbmTimeSlice*) fairManager->GetObject("TimeSlice.");
if (fTimeSlice == NULL) LOG(fatal) << GetName() << ": No time slice branch in the tree!";
} //? time-slice mode
@@ -271,11 +252,11 @@ InitStatus CbmL1::Init()
LOG(info) << GetName() << ": running in event mode.";
- listStsClusters = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsCluster"));
- listStsHitMatch = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHitMatch"));
- listStsClusterMatch = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsClusterMatch"));
+ listStsClusters = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("StsCluster"));
+ listStsHitMatch = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("StsHitMatch"));
+ listStsClusterMatch = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("StsClusterMatch"));
- listStsHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHit"));
+ listStsHits = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("StsHit"));
if (!fUseMUCH) {
fMuchPixelHits = 0;
@@ -288,18 +269,16 @@ InitStatus CbmL1::Init()
listMuchHitMatches = 0;
}
else {
- fMuchPixelHits = (TClonesArray*) fManger->GetObject("MuchPixelHit");
+ fMuchPixelHits = (TClonesArray*) fairManager->GetObject("MuchPixelHit");
}
if (!fUseTRD) {
fTrdPoints = 0;
fTrdHitMatches = 0;
- fTrdPoints = 0;
listTrdHits = 0;
}
else {
-
- listTrdHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("TrdHit"));
+ listTrdHits = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("TrdHit"));
}
if (!fUseTOF) {
@@ -308,42 +287,43 @@ InitStatus CbmL1::Init()
fTofHits = 0;
}
else {
- fTofHits = (TClonesArray*) fManger->GetObject("TofHit");
+ fTofHits = (TClonesArray*) fairManager->GetObject("TofHit");
}
if (fPerformance) {
- CbmMCDataManager* mcManager = (CbmMCDataManager*) fManger->GetObject("MCDataManager");
+ CbmMCDataManager* mcManager = (CbmMCDataManager*) fairManager->GetObject("MCDataManager");
if (NULL == mcManager) LOG(fatal) << GetName() << ": No CbmMCDataManager!";
- fStsPoints = mcManager->InitBranch("StsPoint");
-
fMcEventHeader = mcManager->GetObject("MCEventHeader.");
fMCTracks = mcManager->InitBranch("MCTrack");
- if (NULL == fStsPoints) LOG(fatal) << GetName() << ": No StsPoint data!";
if (NULL == fMCTracks) LOG(fatal) << GetName() << ": No MCTrack data!";
if (NULL == fMcEventHeader) LOG(fatal) << GetName() << ": No MC event header data!";
if (!fLegacyEventMode) {
- fEventList = (CbmMCEventList*) fManger->GetObject("MCEventList.");
+ fEventList = (CbmMCEventList*) fairManager->GetObject("MCEventList.");
if (NULL == fEventList) LOG(fatal) << GetName() << ": No MCEventList data!";
}
if (fUseMVD) {
fMvdPoints = mcManager->InitBranch("MvdPoint");
- listMvdDigiMatches = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdDigiMatch"));
- listMvdHitMatches = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHitMatch"));
+ listMvdDigiMatches = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("MvdDigiMatch"));
+ listMvdHitMatches = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("MvdHitMatch"));
if (!listMvdHitMatches) { LOG(error) << "No listMvdHitMatches provided, performance is not done correctly"; }
}
+ if (fUseSTS) {
+ fStsPoints = mcManager->InitBranch("StsPoint");
+ if (NULL == fStsPoints) LOG(fatal) << GetName() << ": No StsPoint data!";
+ }
+
if (!fUseTRD) {
fTrdPoints = 0;
fTrdHitMatches = 0;
- fTrdPoints = 0;
}
else {
- fTrdHitMatches = (TClonesArray*) fManger->GetObject("TrdHitMatch");
+ fTrdHitMatches = (TClonesArray*) fairManager->GetObject("TrdHitMatch");
fTrdPoints = mcManager->InitBranch("TrdPoint");
}
@@ -353,11 +333,11 @@ InitStatus CbmL1::Init()
}
else {
- fDigisMuch = (TClonesArray*) fManger->GetObject("MuchDigi");
- fDigiMatchesMuch = (TClonesArray*) fManger->GetObject("MuchDigiMatch");
- fClustersMuch = (TClonesArray*) fManger->GetObject("MuchCluster");
+ fDigisMuch = (TClonesArray*) fairManager->GetObject("MuchDigi");
+ fDigiMatchesMuch = (TClonesArray*) fairManager->GetObject("MuchDigiMatch");
+ fClustersMuch = (TClonesArray*) fairManager->GetObject("MuchCluster");
fMuchPoints = mcManager->InitBranch("MuchPoint");
- listMuchHitMatches = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MuchPixelHitMatch"));
+ listMuchHitMatches = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("MuchPixelHitMatch"));
}
if (!fUseTOF) {
@@ -366,14 +346,14 @@ InitStatus CbmL1::Init()
}
else {
fTofPoints = mcManager->InitBranch("TofPoint");
- fTofHitDigiMatches = static_cast<TClonesArray*>(fManger->GetObject("TofHitMatch"));
+ fTofHitDigiMatches = static_cast<TClonesArray*>(fairManager->GetObject("TofHitMatch"));
}
}
else {
}
if (!fUseMVD) { listMvdHits = 0; }
else {
- listMvdHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHit"));
+ listMvdHits = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("MvdHit"));
}
NMvdStationsGeom = 0;
@@ -408,7 +388,6 @@ InitStatus CbmL1::Init()
fpInitManager->InitTargetField(2.5);
-
/**************************************
** **
** STATIONS GEOMETRY INITIALIZATION **
@@ -420,8 +399,10 @@ InitStatus CbmL1::Init()
** Active tracking detector subsystems selection **
***************************************************/
- fActiveTrackingDetectorIDs.insert(L1DetectorID::kSts);
+ fActiveTrackingDetectorIDs.clear();
+
if (fUseMVD) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kMvd); }
+ if (fUseSTS) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kSts); }
if (fUseMUCH) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kMuch); }
if (fUseTRD) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kTrd); }
if (fUseTOF) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kTof); }
@@ -468,8 +449,8 @@ InitStatus CbmL1::Init()
}
}
NTrdStationsGeom = layerCounter;
-
//if ((fTrackingMode == L1Algo::TrackingMode::kMcbm) && (fMissingHits)) { NTrdStationsGeom = NTrdStationsGeom - 1; }
+ LOG(info) << " NTrdStations " << NTrdStationsGeom;
}
/*** ToF ***/
@@ -513,8 +494,16 @@ InitStatus CbmL1::Init()
if (!stsSetup->IsSensorParsInit()) { stsSetup->SetSensorParameters(fStsParSetSensor); }
if (!stsSetup->IsSensorCondInit()) { stsSetup->SetSensorConditions(fStsParSetSensorCond); }
- NMvdStationsGeom = (fUseMVD) ? CbmKF::Instance()->vMvdMaterial.size() : 0;
- NStsStationsGeom = CbmStsSetup::Instance()->GetNofStations();
+ NMvdStationsGeom = 0;
+ if (fUseMVD) {
+ CbmMvdDetector* mvdDetector = CbmMvdDetector::Instance();
+ if (mvdDetector) {
+ CbmMvdStationPar* mvdStationPar = mvdDetector->GetParameterFile();
+ assert(mvdStationPar);
+ NMvdStationsGeom = mvdStationPar->GetStationCount();
+ }
+ }
+ NStsStationsGeom = (fUseSTS) ? CbmStsSetup::Instance()->GetNofStations() : 0;
NStationGeom = NMvdStationsGeom + NStsStationsGeom + NMuchStationsGeom + NTrdStationsGeom + NTOFStationGeom;
// Provide crosscheck number of stations for the fpInitManagera
@@ -592,7 +581,7 @@ InitStatus CbmL1::Init()
// MVD // TODO: to be exchanged with specific flags (timeInfo, fieldInfo etc.) (S.Zh.)
stationInfo.SetTimeInfo(0);
stationInfo.SetTimeResolution(1000.);
- stationInfo.SetFieldStatus(fTrackingMode == L1Algo::TrackingMode::kMcbm ? 0 : 1);
+ stationInfo.SetFieldStatus(L1Algo::TrackingMode::kMcbm == fTrackingMode ? 0 : 1);
stationInfo.SetZ(t.z);
auto thickness = t.dz;
auto radLength = t.RadLength;
@@ -619,7 +608,7 @@ InitStatus CbmL1::Init()
stationInfo.SetStationType(0); // STS
stationInfo.SetTimeInfo(1);
stationInfo.SetTimeResolution(5.);
- stationInfo.SetFieldStatus(fTrackingMode == L1Algo::TrackingMode::kMcbm ? 0 : 1);
+ stationInfo.SetFieldStatus(L1Algo::TrackingMode::kMcbm == fTrackingMode ? 0 : 1);
// Setup station geometry and material
stationInfo.SetZ(cbmSts->GetZ());
double stsXmax = cbmSts->GetXmax();
@@ -697,9 +686,14 @@ InitStatus CbmL1::Init()
fscal trdBackPhi = TMath::Pi() / 2.;
fscal trdFrontSigma = 0.15;
fscal trdBackSigma = 0.15;
+ if (L1Algo::TrackingMode::kGlobal == fTrackingMode) { //SGtrd2D!!
+ trdFrontSigma = 1.1;
+ trdBackSigma = 1.1;
+ stationInfo.SetTimeResolution(1.e10);
+ }
stationInfo.SetFrontBackStripsGeometry(trdFrontPhi, trdFrontSigma, trdBackPhi, trdBackSigma);
stationInfo.SetTrackingStatus(target.z < stationInfo.GetZdouble() ? true : false);
- if (iSt == 1 && fTrackingMode == L1Algo::TrackingMode::kMcbm && fMissingHits) {
+ if (iSt == 1 && L1Algo::TrackingMode::kMcbm == fTrackingMode && fMissingHits) {
stationInfo.SetTrackingStatus(false);
}
fpInitManager->AddStation(stationInfo);
@@ -882,7 +876,7 @@ InitStatus CbmL1::Init()
trackingIterAllSecE.SetTargetPosSigmaXY(10, 10);
// Select default track finder
- if (fTrackingMode == L1Algo::TrackingMode::kMcbm) {
+ if (L1Algo::TrackingMode::kMcbm == fTrackingMode) {
trackingIterAllPrim.SetMaxInvMom(1. / 0.1);
trackingIterAllPrimE.SetMaxInvMom(1. / 0.1);
trackingIterAllSecE.SetMaxInvMom(1. / 0.1);
@@ -936,9 +930,11 @@ void CbmL1::Reconstruct(CbmEvent* event)
{
static int nevent = 0;
vFileEvent.clear();
+ int nStsHits = 0;
+ if (fUseSTS && listStsHits) { nStsHits = listStsHits->GetEntriesFast(); }
L1Vector<std::pair<double, int>> SortStsHits("CbmL1::SortStsHits");
- SortStsHits.reserve(listStsHits->GetEntriesFast());
+ SortStsHits.reserve(nStsHits);
float start_t = 10000000000;
@@ -950,7 +946,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
int FstHitinTs = 0; // 1st hit index in TS
/// sort input hits by time
- for (Int_t j = 0; j < listStsHits->GetEntriesFast(); j++) {
+ for (Int_t j = 0; j < nStsHits; j++) {
CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(j));
double t = sh->GetTime();
if (t < start_t) start_t = t;
@@ -993,11 +989,11 @@ void CbmL1::Reconstruct(CbmEvent* event)
{
int nHits = 0;
int nSta = 1;
- if (listMvdHits) {
+ if (fUseMVD && listMvdHits) {
nHits += listMvdHits->GetEntriesFast();
nSta += NMvdStations;
}
- if (listStsHits) {
+ if (fUseSTS && listStsHits) {
nHits += listStsHits->GetEntriesFast();
nSta += NStsStations;
}
@@ -1128,7 +1124,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
algo->KFTrackFitter_simple();
}
else {
- if (fTrackingMode == L1Algo::TrackingMode::kGlobal || fTrackingMode == L1Algo::TrackingMode::kMcbm) {
+ if (L1Algo::TrackingMode::kGlobal == fTrackingMode || L1Algo::TrackingMode::kMcbm == fTrackingMode) {
algo->L1KFTrackFitterMuch();
}
else {
@@ -1227,8 +1223,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
vRTracksCur.push_back(t);
}
-
- for (int i = 0; i < listStsHits->GetEntriesFast(); i++) {
+ for (int i = 0; i < nStsHits; i++) {
CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(StsIndex[i]));
float time = sh->GetTime();
@@ -1237,7 +1232,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
FstHitinTs = i;
break;
}
- };
+ }
if (!fLegacyEventMode) TsStart = TsStart_new; ///Set new TS strat to earliest discarted track
@@ -1338,7 +1333,6 @@ void CbmL1::Finish()
std::string dir = p.parent_path().string();
if (dir.empty()) dir = ".";
std::string histoOutName = dir + "/L1_histo_" + p.filename().string();
-
LOG(info) << "\033[31;1mHistograms will be saved to: \033[0m" << histoOutName;
TFile* outfile = new TFile(histoOutName.c_str(), "RECREATE");
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 2bd08d0348d624fd5b3c02ce0a566049f497e253..552fabc48ea75fb20c6c5123d4a498a801197729 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -71,6 +71,7 @@ class CbmTofDigiBdfPar;
class CbmTrdParSetDigi;
class CbmTrdParModDigi;
class CbmTofDigiPar;
+class TProfile2D;
class CbmL1HitStore {
public:
@@ -183,13 +184,15 @@ public:
}
/// Correction function for the material budget map
+ /// It fills bins with no statistics
template<L1DetectorID detID>
void ApplyCorrectionToMaterialMap(L1Material& material, const L1MaterialInfo& homogenious)
{
- float hole = 0.;
- if constexpr (detID == L1DetectorID::kMuch || detID == L1DetectorID::kTrd) { hole = 0.15f; }
- else if constexpr (detID == L1DetectorID::kTof) {
- hole = 0.0015f;
+ // TODO: unify the correction function for all detectors
+ float minVal = 0.;
+ if constexpr (detID == L1DetectorID::kMuch) { minVal = 0.15f; }
+ else if constexpr (detID == L1DetectorID::kTof || detID == L1DetectorID::kTrd) {
+ minVal = 0.0015f;
}
// A bit ugly solution, but so can we avoid dependency on input maps file
@@ -197,7 +200,7 @@ public:
if constexpr (detID != L1DetectorID::kSts) { keepRow.resize(material.GetNbins()); }
for (int iBinX = 0; iBinX < material.GetNbins(); ++iBinX) {
- if constexpr (detID == L1DetectorID::kTof) { hole = 0.0015f; }
+ if constexpr (detID == L1DetectorID::kTof) { minVal = 0.0015f; }
if constexpr (detID != L1DetectorID::kSts) {
for (int iBinY = 0; iBinY < material.GetNbins(); ++iBinY) {
keepRow[iBinY] = material.GetRadThick(iBinX, iBinY);
@@ -222,12 +225,18 @@ public:
}
else if constexpr (detID == L1DetectorID::kMuch || detID == L1DetectorID::kTrd || detID == L1DetectorID::kTof) {
// Correction for holes in the material map
- if (iBinY > 0 && iBinY < material.GetNbins() - 1) {
- material.SetRadThick(iBinX, iBinY, TMath::Min(keepRow[iBinY - 1], keepRow[iBinY + 1]));
+ if (L1Algo::TrackingMode::kGlobal != fTrackingMode) {
+ if ((iBinY > 0) && (iBinY < material.GetNbins() - 1)) {
+ material.SetRadThick(iBinX, iBinY, TMath::Min(keepRow[iBinY - 1], keepRow[iBinY + 1]));
+ }
+ }
+ float val = material.GetRadThick(iBinX, iBinY);
+ if (val > 0.0015) { // remember last non-zero value
+ minVal = val;
+ }
+ else { // empty bin with no statistics, fill it with the neoighbours value
+ material.SetRadThick(iBinX, iBinY, minVal);
}
- // Correction for hardcoded values
- if (material.GetRadThick(iBinX, iBinY) > 0.0015) { hole = material.GetRadThick(iBinX, iBinY); }
- if (material.GetRadThick(iBinX, iBinY) < 0.0015) { material.SetRadThick(iBinX, iBinY, hole); }
}
}
}
@@ -235,7 +244,7 @@ public:
/// Utility to map the L1DetectorID items into detector names
- constexpr const char* GetDetectorName(L1DetectorID detectorID)
+ static constexpr const char* GetDetectorName(L1DetectorID detectorID)
{
switch (detectorID) {
case L1DetectorID::kMvd: return "MVD";
@@ -379,6 +388,8 @@ public:
L1Vector<CbmL1HitStore> vHitStore {"CbmL1::vHitStore"}; // diff hit information
private:
+ double LoadMaterialMap(L1Material& mat, const TProfile2D* prof);
+
static CbmL1* fInstance;
L1InitManager* fpInitManager {nullptr}; ///< Pointer to L1InitManager object of L1 algorithm core
@@ -424,6 +435,7 @@ private:
Int_t fTofUseMcHit {-1}; // if Tof data should be processed
Bool_t fUseMVD {false}; // if Mvd data should be processed
+ Bool_t fUseSTS {false}; // if Mvd data should be processed
Bool_t fUseMUCH {false}; // if Much data should be processed
Bool_t fUseTRD {false}; // if Trd data should be processed
Bool_t fUseTOF {false}; // if Tof data should be processed
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index c66d17f13173739587c16ef8ccbc58b98820feb3..adb4636aca070bf5f16d1e574fa36940ad69edf7 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -193,9 +193,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
{ // reserve enough space for hits
int nHitsTotal = 0;
- if (listMvdHits) nHitsTotal += listMvdHits->GetEntriesFast();
+ if (fUseMVD && listMvdHits) nHitsTotal += listMvdHits->GetEntriesFast();
Int_t nEntSts = 0;
- if (listStsHits) {
+ if (fUseSTS && listStsHits) {
if (event) { nEntSts = event->GetNofData(ECbmDataType::kStsHit); }
else {
nEntSts = listStsHits->GetEntriesFast();
@@ -203,9 +203,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
if (nEntSts < 0) nEntSts = 0; // GetNofData() can return -1;
}
nHitsTotal += nEntSts;
- if (fMuchPixelHits) nHitsTotal += fMuchPixelHits->GetEntriesFast();
- if (listTrdHits) nHitsTotal += listTrdHits->GetEntriesFast();
- if (fTofHits) nHitsTotal += fTofHits->GetEntriesFast();
+ if (fUseMUCH && fMuchPixelHits) nHitsTotal += fMuchPixelHits->GetEntriesFast();
+ if (fUseTRD && listTrdHits) nHitsTotal += listTrdHits->GetEntriesFast();
+ if (fUseTOF && fTofHits) nHitsTotal += fTofHits->GetEntriesFast();
tmpHits.reserve(nHitsTotal);
}
@@ -258,8 +258,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
Int_t iFile = set_it->first;
Int_t iEvent = set_it->second;
- // TODO: Probably to put this code in L1Algo interfaces (S.Zharko)
- if (fMvdPoints) {
+ if (fUseMVD && fMvdPoints) {
Int_t nMvdPointsInEvent = fMvdPoints->Size(iFile, iEvent);
double maxDeviation = 0;
for (Int_t iMC = 0; iMC < nMvdPointsInEvent; iMC++) {
@@ -298,7 +297,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // fMvdPoints
firstStsPoint = vMCPoints.size();
- if (fStsPoints) {
+ if (fUseSTS && fStsPoints) {
Int_t nMC = fStsPoints->Size(iFile, iEvent);
double maxDeviation = 0;
for (Int_t iMC = 0; iMC < nMC; iMC++) {
@@ -337,7 +336,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // fStsPoints
firstMuchPoint = vMCPoints.size();
- if (fMuchPoints) {
+ if (fUseMUCH && fMuchPoints) {
Int_t nMC = fMuchPoints->Size(iFile, iEvent);
for (Int_t iMC = 0; iMC < nMC; iMC++) {
CbmL1MCPoint MC;
@@ -365,7 +364,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // fMuchPoints
firstTrdPoint = vMCPoints.size();
- if (fTrdPoints) {
+
+ if (fUseTRD && fTrdPoints) {
for (Int_t iMC = 0; iMC < fTrdPoints->Size(iFile, iEvent); iMC++) {
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 3)) {
@@ -394,7 +394,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
firstTofPoint = vMCPoints.size();
- if (fTofPoints) {
+ if (fUseTOF && fTofPoints) {
vector<float> TofPointToTrackdZ[NTOFStation];
@@ -506,8 +506,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
//
// get MVD hits
- //
- if (listMvdHits) {
+
+ if (fUseMVD && listMvdHits) {
int firstDetStrip = NStrips;
@@ -568,7 +568,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // if listMvdHits
if (fVerbose >= 10) cout << "ReadEvent: mvd hits are gotten." << endl;
- if ((2 == fStsUseMcHit)) { // SG!! create TRD hits from TRD points
+ if (fUseSTS && (2 == fStsUseMcHit)) { // create hits from points
int firstDetStrip = NStrips;
for (int ip = firstStsPoint; ip < firstStsPoint + nStsPoints; ip++) {
@@ -591,7 +591,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
//
// Get STS hits
//
- if (listStsHits && (2 != fStsUseMcHit)) {
+ if (fUseSTS && listStsHits && (2 != fStsUseMcHit)) {
Int_t nEntSts = (event ? event->GetNofData(ECbmDataType::kStsHit) : listStsHits->GetEntriesFast());
@@ -729,7 +729,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
//
// Get MuCh hits
//
- if ((2 == fMuchUseMcHit) && fUseMUCH) { // SG!! create TRD hits from TRD points
+ if ((2 == fMuchUseMcHit) && fUseMUCH) { // create hits from points
int firstDetStrip = NStrips;
for (int ip = firstMuchPoint; ip < firstMuchPoint + nMuchPoints; ip++) {
@@ -836,7 +836,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // for j
} // if listMvdHits
- if ((2 == fTrdUseMcHit) && fTrdHitMatches && fUseTRD) { // SG!! create TRD hits from TRD points
+ if ((2 == fTrdUseMcHit) && fTrdHitMatches && fUseTRD) { // create hits from points
int firstDetStrip = NStrips;
for (int ip = firstTrdPoint; ip < firstTrdPoint + nTrdPoints; ip++) {
@@ -861,12 +861,15 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
int firstDetStrip = NStrips;
vector<bool> mcUsed(nTrdPoints, 0);
- for (int j = 0; j < listTrdHits->GetEntriesFast(); j++) {
+ for (int iHit = 0; iHit < listTrdHits->GetEntriesFast(); iHit++) {
TmpHit th;
- CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(listTrdHits->At(j));
- th.ExtIndex = j;
- th.Det = 3;
+ CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(listTrdHits->At(iHit));
+ // SGtrd2d!! skip TRD-1D hit
+ if ((L1Algo::TrackingMode::kGlobal == fTrackingMode) && (int) mh->GetClassType() != 1) { continue; }
+
+ th.ExtIndex = iHit;
+ th.Det = 3;
th.id = tmpHits.size();
@@ -884,7 +887,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.dt = mh->GetTimeError();
// th.iSector = 0;
- th.iStripF = firstDetStrip + j;
+ th.iStripF = firstDetStrip + iHit;
th.iStripB = th.iStripF; //TrdHitsOnStationIndex[num+1][k];
if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
@@ -908,73 +911,71 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
- th.iMC = -1;
- int iMC = -1;
+ th.iMC = -1;
if (fPerformance && fTrdHitMatches) {
- CbmMatch* trdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(j));
+ CbmMatch* trdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(iHit));
+
+ if (trdHitMatch->GetNofLinks() > 0) {
- if (trdHitMatch->GetNofLinks() > 0)
- if (trdHitMatch->GetLink(0).GetIndex() < nTrdPoints) {
- iMC = trdHitMatch->GetLink(0).GetIndex();
- th.iMC = iMC + nMvdPoints + nStsPoints + nMuchPoints;
- // th.track = vMCPoints[th.iMC].ID;
+ int iMC = trdHitMatch->GetLink(0).GetIndex();
- if ((1 == fTrdUseMcHit) && (th.iMC > -1)) //SG!!! replace hits by MC points
+ assert(iMC >= 0 && iMC < nTrdPoints);
- th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetStations()[th.iStation]);
+ th.iMC = iMC + nMvdPoints + nStsPoints + nMuchPoints;
+ th.track = vMCPoints[th.iMC].ID;
- if (L1Algo::TrackingMode::kGlobal == fTrackingMode) { //SG!!! replace hits by MC points
+ if (1 == fTrdUseMcHit) { // replace hits by MC points
+ if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
+ // SGtrd2d!! don't smear
+ th.dx = 0.;
+ th.dy = 0.;
+ th.dt = 0.;
+ }
+ th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetStations()[th.iStation]);
+ /*
CbmTrdPoint* pt =
(CbmTrdPoint*) fTrdPoints->Get(trdHitMatch->GetLink(0).GetFile(), trdHitMatch->GetLink(0).GetEntry(),
trdHitMatch->GetLink(0).GetIndex());
-
- th.dx = 0.1;
- th.du = 0.1;
- th.dy = 0.1;
- th.dv = 0.1;
-
- th.x = 0.5 * (pt->GetXOut() + pt->GetXIn()); // + gRandom->Gaus(0, th.dx);
- th.y = 0.5 * (pt->GetYOut() + pt->GetYIn()); // + gRandom->Gaus(0, th.dy);
- th.z = 0.5 * (pt->GetZOut() + pt->GetZIn());
- th.time = pt->GetTime();
- th.dt = 10000.;
-
- th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
- th.dx = 1.;
- th.du = 1.;
- th.dy = 1.;
- th.dv = 1.;
-
- if (th.iMC < 0) continue;
- if (mcUsed[iMC]) continue;
- if (0) {
- int mcTrack = -1;
- float mcP = 0;
- if (th.iMC >= 0) {
- mcTrack = vMCPoints[th.iMC].ID;
- if (mcTrack >= 0) { mcP = vMCTracks[mcTrack].p; }
- }
- if (mcP < 1.) continue;
+
+ cout << " dx " << th.x - 0.5 * (pt->GetXOut() + pt->GetXIn()) << " dy "
+ << th.y - 0.5 * (pt->GetYOut() + pt->GetYIn()) << " dz "
+ << th.z - 0.5 * (pt->GetZOut() + pt->GetZIn()) << " z " << th.z << endl;
+ */
+ if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
+ // SGtrd2d!! enlarge errors
+ th.dx = 1.1;
+ th.du = 1.1;
+ th.dy = 1.1;
+ th.dv = 1.1;
+ th.dt = 10000.;
+ }
+ if (th.iMC < 0) continue;
+ if (mcUsed[iMC]) continue;
+ if (0) { // select specific tracks only
+ int mcTrack = -1;
+ float mcP = 0;
+ if (th.iMC >= 0) {
+ mcTrack = vMCPoints[th.iMC].ID;
+ if (mcTrack >= 0) { mcP = vMCTracks[mcTrack].p; }
}
- mcUsed[iMC] = 1;
+ if (mcP < 1.) continue;
}
+ mcUsed[iMC] = 1;
}
+ }
}
-
tmpHits.push_back(th);
nTrdHits++;
-
- } // for j
+ } // for listTrdHits
} // if listTrdHits
//
// Get ToF hits
//
- if ((2 == fTofUseMcHit) && fUseTOF) { // SG!! create TRD hits from TRD points
+ if ((2 == fTofUseMcHit) && fUseTOF) { // create hits from points
int firstDetStrip = NStrips;
for (int ip = firstTofPoint; ip < firstTofPoint + nTofPoints; ip++) {
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 1b1ec5a1c5fa467a9b461a82d28ff7574ac85ec5..6e96d7217b457d863514e8c2631e09700d0e697f 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -240,7 +240,7 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
T.C40 = T.C41 = T.C42 = T.C43 = 0;
T.C50 = T.C51 = T.C52 = T.C53 = T.C54 = 0;
T.C22 = T.C33 = fMaxSlopePV * fMaxSlopePV / 9.;
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) T.C22 = T.C33 = 10;
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) T.C22 = T.C33 = 10;
T.C44 = fMaxInvMom / 3. * fMaxInvMom / 3.;
T.C55 = timeEr * timeEr;
@@ -262,10 +262,8 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
}
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm)
- // add the target
- {
-
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
+ // add the target
if (istal < fNfieldStations) {
fvec eX, eY, J04, J14;
fvec dz;
@@ -315,15 +313,15 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
T.C00 = TargetXYInfo.C00;
T.C10 = TargetXYInfo.C10;
T.C11 = TargetXYInfo.C11;
-
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) // extrapolate to left hit
- {
+ // extrapolate to left hit
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
if (istal <= fNfieldStations) L1Extrapolate0(T, zl, fld0);
else
L1ExtrapolateLine(T, zl);
}
- else
+ else {
L1Extrapolate0(T, zl, fld0);
+ }
for (int ista = 0; ista <= istal - 1; ista++) {
if constexpr (L1Parameters::kIfUseRadLengthTable) {
@@ -340,33 +338,39 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
if (fUseHitErrors) info.sigma2 = du1 * du1;
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
if (istal < fNfieldStations) L1Filter(T, info, u);
- else
+ else {
L1FilterNoField(T, info, u);
+ }
}
- else
+ else {
L1Filter(T, info, u);
-
+ }
info = stal.backInfo;
if (fUseHitErrors) info.sigma2 = dv1 * dv1;
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
if (istal < fNfieldStations) L1Filter(T, info, v);
- else
+ else {
L1FilterNoField(T, info, v);
+ }
}
- else
+ else {
L1Filter(T, info, v);
+ }
#endif // BEGIN_FROM_TARGET
if constexpr (L1Parameters::kIfUseRadLengthTable) {
- if (fTrackingMode != kMcbm) fit.L1AddMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), fMaxInvMom, 1);
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm)
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
fit.L1AddThickMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), fMaxInvMom, 1, stal.materialInfo.thick, 1);
+ }
+ else {
+ fit.L1AddMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), fMaxInvMom, 1);
+ }
}
else {
fit.L1AddMaterial(T, stal.materialInfo, fMaxInvMom, 1);
@@ -378,16 +382,16 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
fvec dz = zstam - zl;
L1ExtrapolateTime(T, dz, stam.timeInfo);
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm)
// extrapolate to left hit
- {
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
if (istam < fNfieldStations) L1Extrapolate0(T, zstam, fld0);
- else
+ else {
L1ExtrapolateLine(T, zstam); // TODO: fld1 doesn't work!
+ }
}
- else
+ else {
L1Extrapolate0(T, zstam, fld0); // TODO: fld1 doesn't work!
-
+ }
} // i1_V
}
@@ -430,9 +434,9 @@ inline void L1Algo::f20(
L1HitIndex_t imh = 0;
- int irm1 = -1;
+ int irm1 = -1;
while (true) {
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
irm1++;
if ((L1HitIndex_t) irm1
>= (StsHitsUnusedStopIndex[&stam - fStations.begin()] - StsHitsUnusedStartIndex[&stam - fStations.begin()]))
@@ -442,8 +446,8 @@ inline void L1Algo::f20(
else {
if (!areaTime.GetNext(imh)) break;
}
- const L1HitPoint& hitm = vStsHits_m[imh];
+ const L1HitPoint& hitm = vStsHits_m[imh];
// check y-boundaries
if ((stam.timeInfo) && (stal.timeInfo))
@@ -663,35 +667,38 @@ inline void L1Algo::f30( // input
if (fUseHitErrors) info.sigma2 = du2 * du2;
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
if (istam < fNfieldStations) L1Filter(T2, info, u_front_2);
- else
+ else {
L1FilterNoField(T2, info, u_front_2);
+ }
}
- else
+ else {
L1Filter(T2, info, u_front_2);
-
+ }
info = stam.backInfo;
if (fUseHitErrors) info.sigma2 = dv2 * dv2;
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
if (istam < fNfieldStations) L1Filter(T2, info, u_back_2);
- else
+ else {
L1FilterNoField(T2, info, u_back_2);
+ }
}
- else
+ else {
L1Filter(T2, info, u_back_2);
-
+ }
FilterTime(T2, timeM, timeMEr, stam.timeInfo);
if constexpr (L1Parameters::kIfUseRadLengthTable) {
- if (fTrackingMode != kMcbm && fTrackingMode != kGlobal)
- fit.L1AddMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp, 1);
-
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm)
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
fit.L1AddThickMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), fMaxInvMom, 1, stam.materialInfo.thick,
1);
+ }
+ else {
+ fit.L1AddMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp, 1);
+ }
}
else {
fit.L1AddMaterial(T2, stam.materialInfo, T2.qp, 1);
@@ -705,19 +712,21 @@ inline void L1Algo::f30( // input
// extrapolate to the right hit station
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm)
- // extrapolate to the right hit station
- {
- if (istar <= fNfieldStations) L1Extrapolate(T2, star.z, T2.qp, f2); // Full extrapolation in the magnetic field
- else
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
+ if (istar <= fNfieldStations) {
+ L1Extrapolate(T2, star.z, T2.qp, f2); // Full extrapolation in the magnetic field
+ }
+ else {
L1ExtrapolateLine(T2, star.z); // Extrapolation with line ()
+ }
}
- else
+ else {
L1Extrapolate(T2, star.z, T2.qp, f2);
+ }
// ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----
for (Tindex i2_4 = 0; i2_4 < n2_4; ++i2_4) {
- if (fTrackingMode == kSts && (T2.C44[i2_4] < 0)) { continue; }
+ if (kSts == fTrackingMode && (T2.C44[i2_4] < 0)) { continue; }
if (T2.C00[i2_4] < 0 || T2.C11[i2_4] < 0 || T2.C22[i2_4] < 0 || T2.C33[i2_4] < 0 || T2.C55[i2_4] < 0) continue;
if (fabs(T2.tx[i2_4]) > fMaxSlope) continue;
@@ -739,9 +748,9 @@ inline void L1Algo::f30( // input
L1HitIndex_t irh = 0;
L1HitIndex_t Ntriplets = 0;
- int irh1 = -1;
+ int irh1 = -1;
while (true) {
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
irh1++;
if ((L1HitIndex_t) irh1 >= (StsHitsUnusedStopIndex[istar] - StsHitsUnusedStartIndex[istar])) break;
irh = irh1;
@@ -929,29 +938,26 @@ inline void L1Algo::f31( // input
if (fUseHitErrors) info.sigma2 = du_[i3_V] * du_[i3_V];
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
- if ((&star - fStations.begin()) < fNfieldStations) L1Filter(T_3[i3_V], info, u_front_[i3_V]);
- else {
- L1FilterNoField(T_3[i3_V], info, u_front_[i3_V]);
- }
- }
- else
- L1Filter(T_3[i3_V], info, u_front_[i3_V]);
+ bool noField =
+ (kGlobal == fTrackingMode || kMcbm == fTrackingMode) && (&star - fStations.begin() >= fNfieldStations);
+ if (noField) { L1FilterNoField(T_3[i3_V], info, u_front_[i3_V]); }
+ else {
+ L1Filter(T_3[i3_V], info, u_front_[i3_V]);
+ }
info = star.backInfo;
if (fUseHitErrors) info.sigma2 = dv_[i3_V] * dv_[i3_V];
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
- if ((&star - fStations.begin()) < fNfieldStations) L1Filter(T_3[i3_V], info, u_back_[i3_V]);
- else
- L1FilterNoField(T_3[i3_V], info, u_back_[i3_V]);
- }
- else
+ if (noField) { L1FilterNoField(T_3[i3_V], info, u_back_[i3_V]); }
+ else {
L1Filter(T_3[i3_V], info, u_back_[i3_V]);
+ }
- if (fTrackingMode != kMcbm) FilterTime(T_3[i3_V], timeR[i3_V], timeER[i3_V], star.timeInfo);
+ if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) {
+ FilterTime(T_3[i3_V], timeR[i3_V], timeER[i3_V], star.timeInfo);
+ }
// FilterTime(T_3[i3_V], timeR[i3_V], timeER[i3_V]);
}
@@ -1535,8 +1541,8 @@ L1Algo::TripletsStaPort( /// creates triplets: input: @istal - start station nu
L1HitIndex_t* RealIHitM = &((*RealIHitP)[StsHitsUnusedStartIndex[istam]]);
L1HitIndex_t* RealIHitR = &((*RealIHitP)[StsHitsUnusedStartIndex[istar]]);
for (Tindex i = 0; i < n3; ++i) {
- Tindex i_4 = i % 4;
- Tindex i_V = i / 4;
+ Tindex i_4 = i % 4;
+ Tindex i_V = i / 4;
L1HitIndex_t iHits[3] = {RealIHitL[hitsl_3[i]], RealIHitM[hitsm_3[i]], RealIHitR[hitsr_3[i]]};
#ifdef PULLS
if (fL1Eff_triplets->AddOne(iHits)) fL1Pulls->AddOne(T_3[i_V], i_4, iHits[2]);
@@ -1659,8 +1665,8 @@ void L1Algo::CATrackFinder()
static Tindex stat_nTriplets[fNFindIterations] = {0};
static Tindex stat_nLevels[L1Parameters::kMaxNstations - 2][fNFindIterations] = {{0}};
- static Tindex stat_nCalls[fNFindIterations] = {0}; // n calls of CAFindTrack
- static Tindex stat_nTrCandidates[fNFindIterations] = {0};
+ static Tindex stat_nCalls[fNFindIterations] = {0}; // n calls of CAFindTrack
+ static Tindex stat_nTrCandidates[fNFindIterations] = {0};
#endif
/********************************/ /**
@@ -1778,7 +1784,7 @@ void L1Algo::CATrackFinder()
// ---- Loop over Track Finder iterations ----------------------------------------------------------------//
L1ASSERT(0, fNFindIterations == fParameters.CAIterationsContainer().size());
- isec = 0; // TODO: temporary! (S.Zharko)
+ isec = 0; // TODO: temporary! (S.Zharko)
for (const auto& caIteration : fParameters.CAIterationsContainer()) // all finder
{
std::cout << "CA Track Finder Iteration!!" << isec << '\n';
@@ -1799,8 +1805,8 @@ void L1Algo::CATrackFinder()
if (isec != 0) {
L1Vector<L1HitIndex_t>* RealIHitPTmp = RealIHitP;
- RealIHitP = RealIHitPBuf;
- RealIHitPBuf = RealIHitPTmp;
+ RealIHitP = RealIHitPBuf;
+ RealIHitPBuf = RealIHitPTmp;
L1Vector<L1Hit>* vStsHitsUnused_temp = vStsHitsUnused;
vStsHitsUnused = vStsHitsUnused_buf;
@@ -1828,7 +1834,6 @@ void L1Algo::CATrackFinder()
// --- SET PARAMETERS FOR THE ITERATION ---
FIRSTCASTATION = 0;
- if (fTrackingMode == kGlobal) { FIRSTCASTATION = 12; }
// if ( (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) )
// FIRSTCASTATION = 2;
@@ -2214,9 +2219,10 @@ void L1Algo::CATrackFinder()
if (first_trip.GetLevel() == 0)
continue; // ghost suppression // find track with 3 hits only if it was created from a chain of triplets, but not from only one triplet
- if (fTrackingMode != kMcbm)
+ if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) {
if ((ilev == 0) && (GetFStation((*fStripFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]) != 0))
continue; // ghost supression // collect only MAPS tracks-triplets CHECK!!!
+ }
}
if (first_trip.GetLevel() < ilev)
continue; // try only triplets, which can start track with ilev+3 length. w\o it have more ghosts, but efficiency either
@@ -2360,7 +2366,10 @@ void L1Algo::CATrackFinder()
}
}
- if (fTrackingMode == kMcbm) {
+ if (kMcbm == fTrackingMode) {
+ if (tr.NHits <= 3) { check = 0; }
+ }
+ else if (kGlobal == fTrackingMode) {
if (tr.NHits <= 3) { check = 0; }
}
else {
@@ -2369,8 +2378,9 @@ void L1Algo::CATrackFinder()
if (check) {
#ifdef EXTEND_TRACKS
- if (fTrackingMode != kMcbm)
+ if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) {
if (tr.NHits != fNstations) BranchExtender(tr);
+ }
#endif
float sumTime = 0;
@@ -2385,6 +2395,9 @@ void L1Algo::CATrackFinder()
#endif
+ if (kGlobal == fTrackingMode) { //SGtrd2d!!
+ if (tr.fStsHits.size() < 4) continue;
+ }
for (L1Vector<L1HitIndex_t>::iterator phIt = tr.fStsHits.begin(); /// used strips are marked
phIt != tr.fStsHits.end(); ++phIt) {
L1Hit& h = (*vStsHits)[*phIt];
@@ -2467,11 +2480,11 @@ void L1Algo::CATrackFinder()
int NHitsOnStation = 0;
for (int ista = 0; ista < fNstations; ++ista) {
- int start = StsHitsUnusedStartIndex[ista];
- int Nelements = StsHitsUnusedStopIndex[ista] - start;
- L1Hit* staHits = nullptr; // to avoid out-of-range error in ..[start]
+ int start = StsHitsUnusedStartIndex[ista];
+ int Nelements = StsHitsUnusedStopIndex[ista] - start;
+ L1Hit* staHits = nullptr; // to avoid out-of-range error in ..[start]
L1HitIndex_t* staHitIndices = nullptr;
- L1HitPoint* staHitPoints = nullptr;
+ L1HitPoint* staHitPoints = nullptr;
if (Nelements > 0) {
staHits = &((*vStsHitsUnused)[start]);
staHitIndices = &((*RealIHitP)[start]);
@@ -2723,7 +2736,7 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
fscal dqp = fabs(qp1 - qp2);
fscal Cqp = curr_trip->GetCqp();
Cqp += new_trip.GetCqp();
- if (fTrackingMode != kMcbm) {
+ if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) {
if (dqp > fPickNeighbour * Cqp) {
continue; // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
}
@@ -2741,11 +2754,16 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
fscal Cty = curr_trip->GetCty();
Cty += new_trip.GetCty();
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
+ if (kMcbm == fTrackingMode) {
if (dty > 3 * Cty) continue;
if (dtx > 3 * Ctx) continue;
}
+ if (kGlobal == fTrackingMode) {
+ if (dty > 6 * sqrt(Cty)) continue; //SGtrd2d
+ if (dtx > 7 * sqrt(Ctx)) continue;
+ }
+
if (GetFUsed((*fStripFlag)[(*vStsHitsUnused)[new_trip.GetLHit()].f]
| (*fStripFlag)[(*vStsHitsUnused)[new_trip.GetLHit()].b])) { //hits are used
// no used hits allowed -> compare and store track
@@ -2770,17 +2788,29 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
dqp = dqp / Cqp;
- dtx = dtx / Ctx; // TODO: SG: it must be /sqrt(Ctx);
- dty = dty / Cty; // TODO: SG: it must be /sqrt(Cty);
+ if (kGlobal == fTrackingMode) { //SGtrd2d!!!
+ dtx = dtx / sqrt(Ctx);
+ dty = dty / sqrt(Cty);
+ }
+ else {
+ dtx = dtx / Ctx; // TODO: SG: it must be /sqrt(Ctx);
+ dty = dty / Cty; // TODO: SG: it must be /sqrt(Cty);
+ }
- if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
new_chi2 += dtx * dtx;
new_chi2 += dty * dty;
}
else {
new_chi2 += dqp * dqp;
}
- if (new_chi2 > fTrackChi2Cut * new_L) continue; // TODO: SG: it must be ( 2 * new_L )
+
+ if (kGlobal == fTrackingMode) { //SGtrd2d!!!
+ if (new_chi2 > fTrackChi2Cut * 2 * new_L) continue;
+ }
+ else {
+ if (new_chi2 > fTrackChi2Cut * new_L) continue; // TODO: SG: it must be ( 2 * new_L )
+ }
const int new_ista = ista + new_trip.GetMSta() - new_trip.GetLSta();
CAFindTrack(new_ista, best_tr, best_L, best_chi2, &new_trip, new_tr[ista], new_L, new_chi2, min_best_l, new_tr);
diff --git a/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx b/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
index 128364d7730fe792027730fe020f0fa506f5db32..c2a399025f3ca832a3bbd373e04c6e11c5bb5e79 100644
--- a/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
+++ b/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
@@ -63,7 +63,6 @@ void CbmL1StsTrackFinder::Init() {}
// ----- Copy tracks to output array -----------------------------------
Int_t CbmL1StsTrackFinder::CopyL1Tracks(CbmEvent* event)
{
-
CbmL1* L1 = CbmL1::Instance();
if (!L1) return 0;
@@ -91,6 +90,9 @@ Int_t CbmL1StsTrackFinder::CopyL1Tracks(CbmEvent* event)
for (vector<int>::iterator ih = it->StsHits.begin(); ih != it->StsHits.end(); ++ih) {
CbmL1HitStore& h = L1->vHitStore[*ih];
// double zref = L1->algo->vStations[h.iStation].z[0];
+ if (h.Det > 1) { // not MVD or STS hit
+ continue;
+ }
if (h.ExtIndex < 0) {
// CbmMvdHit tmp;
// tmp.SetZ(zref);
diff --git a/reco/detectors/trd/CbmTrdModuleRec2D.cxx b/reco/detectors/trd/CbmTrdModuleRec2D.cxx
index 55f3095709c285cf6f1342a0aad0397704f9b7ad..d6ed10c2d04047d0388982b876f2cced20de15f0 100644
--- a/reco/detectors/trd/CbmTrdModuleRec2D.cxx
+++ b/reco/detectors/trd/CbmTrdModuleRec2D.cxx
@@ -1022,7 +1022,7 @@ Int_t CbmTrdModuleRec2D::LoadDigis(vector<const CbmTrdDigi*>* din, Int_t cid)
//row =
GetPadRowCol(dgR->GetAddressChannel(), colR);
}
- if (colR == colT) {
+ if (colR == colT && dgR != NULL) {
fDigis[cid].push_back(new CbmTrdDigiRec(*dgT, *dgR));
j = din->erase(j);
}