diff --git a/macro/beamtime/common/CMakeLists.txt b/macro/beamtime/common/CMakeLists.txt index 72153848805140d9c87c73810a62af4569d682ad..62107b9d4a4ed6883aaf94a04638fd8f96b1efda 100644 --- a/macro/beamtime/common/CMakeLists.txt +++ b/macro/beamtime/common/CMakeLists.txt @@ -24,3 +24,12 @@ Install(DIRECTORY online Install(DIRECTORY data DESTINATION share/cbmroot/macro/beamtime/common PATTERN "*" EXCLUDE) + +If(DEFINED ENV{RAW_DATA_PATH} ) + # ===== Copy the .rootrc file into the directory from which root is executed + # --- Otherwise the rootalias file is not loaded + file(COPY ${CBMROOT_SOURCE_DIR}/macro/include/.rootrc + DESTINATION ${CBMROOT_BINARY_DIR}/macro/beamtime/common) + +EndIf() # If(DEFINED ENV{RAW_DATA_PATH} ) + diff --git a/macro/beamtime/mcbm2025/mcbm_event_reco_L1.C b/macro/beamtime/mcbm2025/mcbm_event_reco_L1.C new file mode 100644 index 0000000000000000000000000000000000000000..ee0f42f363fb8a34cc9a78275877d580be89d7d9 --- /dev/null +++ b/macro/beamtime/mcbm2025/mcbm_event_reco_L1.C @@ -0,0 +1,738 @@ +/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file mcbm_event_reco_L1.C +/// @brief Offline reconstruction of real data in mCBM (2025) +/// @since 14.02.2025 +/// @author Sergei Zharko <s.zharko@gsi.de> + + +#include <cmath> +#include <cstdio> +#include <string> + +/// FIXME: Disable clang formatting to keep easy parameters overview +/* clang-format off */ +Bool_t mcbm_event_reco_L1(UInt_t uRunId = 3453, + Int_t nTimeslices = 10, + TString sInpDir = "./data/", + TString sOutDir = "./data/", + Int_t iUnpFileIndex = -1, + Bool_t bMVD = kFALSE, + Bool_t bSTS = kTRUE, + Bool_t bTRD = kTRUE, + Bool_t bTRD2d = kTRUE, + Bool_t bRICH = kTRUE, + Bool_t bMUCH = kFALSE, + Bool_t bTOF = kTRUE, + Bool_t bTOFtr = kFALSE, + Bool_t bPSD = kFALSE, + Bool_t bAli = kTRUE, + Bool_t bEvB = kTRUE, + Bool_t bL1 = kTRUE, + Bool_t bQA = kTRUE, + Bool_t bFSD = kFALSE, + TString sInpFile = "", + Bool_t bPV = kFALSE + ) +{ + /// FIXME: Re-enable clang formatting after parameters initial values setting + /* clang-format on */ + + // --- Logger settings ---------------------------------------------------- + TString logLevel = "INFO"; //"INFO"; + TString logVerbosity = "LOW"; //"VERYLOW"; + // ------------------------------------------------------------------------ + + // ----- Environment -------------------------------------------------- + TString myName = "mcbm_event_reco_L1"; // this macro's name for screen output + TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory + // ------------------------------------------------------------------------ + + + // ----- In- and output file names ------------------------------------ + /// Standardized RUN ID + TString sRunId = TString::Format("%04u", uRunId); + + gSystem->Exec("mkdir -p " + sOutDir); + gSystem->Exec("cp $VMCWORKDIR/macro/run/.rootrc ."); + + TString outFile = sOutDir + "/" + sRunId; + /// Initial pattern + TString inFile = sInpDir + sRunId + ".digi"; + + /// Add index of splitting at unpacking level if needed + if (0 <= iUnpFileIndex) { + inFile += TString::Format("_%02u", iUnpFileIndex); + // the input par file is not split during unpacking! + outFile += TString::Format("_%02u", iUnpFileIndex); + } // if ( 0 <= uUnpFileIndex ) + /// Add ROOT file suffix + inFile += ".root"; + TString parFileOut = outFile + ".par.root"; + TString geoFileOut = outFile + ".geo.root"; + outFile += ".rec.root"; + + if ("" != sInpFile) { + inFile = sInpFile; + outFile = inFile; + outFile.ReplaceAll(".digi.root", ""); + parFileOut = outFile + ".par.root"; + geoFileOut = outFile + ".geo.root"; + outFile += ".rec.root"; + } + // --------------------------------------------- + + + // ----- TOF defaults ------------------------ + // Your folder with the Tof Calibration files; + TString cFileId = sRunId + ".5.000"; + TString TofFileFolder = ""; + + // ===> PAL 2022/11/04: overwriten by block around l.510! + Int_t calMode = 93; + Int_t calSel = 1; + Int_t calSm = 2; + Int_t RefSel = 11; + Double_t dDeadtime = 50.; + Int_t iSel2 = -1; + + // Tracking + Int_t iSel = 22002; + Int_t iTrackingSetup = 1; // 2 for checking without beam counter; + Int_t iGenCor = 1; + Double_t dScalFac = 2.5; + Double_t dChi2Lim2 = 4.; + Bool_t bUseSigCalib = kFALSE; + Int_t iCalOpt = 110; // 0=do not call CbmTofCalibrator + Int_t iTrkPar = 0; // 4 for check without beam counter + Double_t dTOffScal = 1.; + // ------------------------------------------------------------------------ + + // ----- TOF Calibration Settings --------------------------------------- + TString cCalId = "490.100.5.0"; + if (uRunId >= 759) cCalId = "759.100.4.0"; + if (uRunId >= 812) cCalId = "831.100.4.0"; + if (uRunId >= 1588) cCalId = "1588.50.6.0"; + if (uRunId >= 2160) cCalId = "2160.50.4.0"; + if (uRunId >= 2352) cCalId = "2391_1"; + if (uRunId >= 2700) cCalId = "2912.1"; + if (uRunId >= 2900) cCalId = "3026_1"; + if (uRunId >= 3399) cCalId = "3310_1"; + + Int_t iCalSet = 30040500; // calibration settings + if (uRunId >= 759) iCalSet = 10020500; + if (uRunId >= 812) iCalSet = 10020500; + if (uRunId >= 1588) iCalSet = 12002002; + if (uRunId >= 2160) iCalSet = 700900500; + if (uRunId >= 2352) iCalSet = 22002500; + if (uRunId >= 2700) iCalSet = 12032500; + + Double_t Tint = 100.; // coincidence time interval + Int_t iTrackMode = 2; // 2 for TofTracker + const Int_t iTofCluMode = 1; + // ------------------------------------------------------------------------ + + // --- Load the geometry setup ---- + // This is currently only required by the TRD (parameters) + cbm::mcbm::ToForceLibLoad dummy; /// Needed to trigger loading of the library as no fct dict in ROOT6 and CLING + TString geoSetupTag = ""; + try { + geoSetupTag = cbm::mcbm::GetSetupFromRunId(uRunId); + } + catch (const std::invalid_argument& e) { + std::cout << "Error in mapping from runID to setup name: " << e.what() << std::endl; + return kFALSE; + } + + TString geoFile = sInpDir + "/" + geoSetupTag + ".geo.root"; + CbmSetup* geoSetup = CbmSetup::Instance(); + geoSetup->LoadSetup(geoSetupTag); + + // You can modify the pre-defined setup by using + geoSetup->SetActive(ECbmModuleId::kMvd, bMVD); + geoSetup->SetActive(ECbmModuleId::kSts, bSTS); + geoSetup->SetActive(ECbmModuleId::kMuch, bMUCH); + geoSetup->SetActive(ECbmModuleId::kRich, bRICH); + geoSetup->SetActive(ECbmModuleId::kTrd, bTRD); + geoSetup->SetActive(ECbmModuleId::kTrd2d, bTRD2d); + geoSetup->SetActive(ECbmModuleId::kTof, bTOF); + geoSetup->SetActive(ECbmModuleId::kPsd, bPSD); + geoSetup->SetActive(ECbmModuleId::kFsd, bFSD); + + + //----- Load Parameters -------------------------------------------------- + TList* parFileList = new TList(); + TofFileFolder = Form("%s/%s", TofFileFolder.Data(), cCalId.Data()); + + // ----- TRD digitisation parameters ------------------------------------- + TString geoTagTrd; + if (geoSetup->IsActive(ECbmModuleId::kTrd)) { + if (geoSetup->GetGeoTag(ECbmModuleId::kTrd, geoTagTrd)) { + TString paramFilesTrd(Form("%s/parameters/trd/trd_%s", srcDir.Data(), geoTagTrd.Data())); + std::vector<TString> paramFilesVecTrd = {"asic", "digi", "gas", "gain"}; + for (auto parIt : paramFilesVecTrd) { + parFileList->Add(new TObjString(Form("%s.%s.par", paramFilesTrd.Data(), parIt.Data()))); + } + } + for (auto parFileVecIt : *parFileList) { + LOG(debug) << Form("TrdParams - %s - added to parameter file list\n", parFileVecIt->GetName()); + } + } + // ----- TOF digitisation parameters ------------------------------------- + TString geoTagTof; + if (geoSetup->IsActive(ECbmModuleId::kTof)) { + geoSetup->GetGeoTag(ECbmModuleId::kTof, geoTagTof); + TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTagTof + ".digibdf.par"); + parFileList->Add(tofBdfFile); + std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl; + + // TFile* fgeo = new TFile(geoFile); + // TGeoManager* geoMan = (TGeoManager*) fgeo->Get("FAIRGeom"); + // if (NULL == geoMan) { + // cout << "<E> FAIRGeom not found in geoFile " << geoFile.Data() << endl; + // return 1; + // } + } + // ------------------------------------------------------------------------ + + // ----- Timer -------------------------------------------------------- + TStopwatch timer; + timer.Start(); + // ------------------------------------------------------------------------ + + + // ----- FairRunAna --------------------------------------------------- + FairRunAna* run = new FairRunAna(); + FairFileSource* inputSource = new FairFileSource(inFile); + run->SetSource(inputSource); + + FairRootFileSink* outputSink = new FairRootFileSink(outFile); + run->SetSink(outputSink); + run->SetGeomFile(geoFile); + + // Define output file for FairMonitor histograms + TString monitorFile{outFile}; + monitorFile.ReplaceAll(".rec.root", ".rec.monitor.root"); + FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile); + // ------------------------------------------------------------------------ + + + // ----- Logger settings ---------------------------------------------- + FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data()); + FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data()); + //FairLogger::GetLogger()->SetLogScreenLevel("DEBUG"); + // ------------------------------------------------------------------------ + + + // ========================================================================= + // === Alignment Correction === + // ========================================================================= + // (Fairsoft Apr21p2 or newer is needed) + if (bAli) { + + TString alignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geoSetupTag + ".root"; + + if (alignmentMatrixFileName.Length() != 0) { + + std::cout << "-I- " << myName << ": Applying alignment for file " << alignmentMatrixFileName << std::endl; + + // Define the basic structure which needs to be filled with information + // This structure is stored in the output file and later passed to the + // FairRoot framework to do the (miss)alignment + std::map<std::string, TGeoHMatrix>* matrices{nullptr}; + + // read matrices from disk + LOG(info) << "Filename: " << alignmentMatrixFileName; + TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ"); + if (misalignmentMatrixRootfile->IsOpen()) { + gDirectory->GetObject("MisalignMatrices", matrices); + misalignmentMatrixRootfile->Close(); + } + else { + LOG(error) << "Could not open file " << alignmentMatrixFileName << "\n Exiting"; + exit(1); + } + + if (matrices) { + run->AddAlignmentMatrices(*matrices); + } + else { + LOG(error) << "Alignment required but no matrices found." + << "\n Exiting"; + exit(1); + } + } + } + // ------------------------------------------------------------------------ + + + // --------------------event builder--------------------------------------- + // ----- EventBuilder Settings---------------- + // mCbm track trigger Tof, Bmon & STS (case 4 of mcbm_unp_event.C) + const Int_t eb_TriggerMinNumberBmon{1}; + const Int_t eb_TriggerMaxNumberBmon{2}; + const Int_t eb_TriggerMinNumberSts{2}; + const Int_t eb_TriggerMinNumberStsLayers{1}; + const Int_t eb_TriggerMinNumberMuch{1}; + const Int_t eb_TriggerMinNumberTof{8}; + const Int_t eb_TriggerMinNumberTofLayers{4}; + const Int_t eb_TriggerMinNumberRich{0}; + + const Int_t eb_TrigWinMinBmon{-50}; + const Int_t eb_TrigWinMaxBmon{50}; + const Int_t eb_TrigWinMinSts{-60}; + const Int_t eb_TrigWinMaxSts{60}; + const Int_t eb_TrigWinMinMuch{-100}; + const Int_t eb_TrigWinMaxMuch{100}; + const Int_t eb_TrigWinMinTrd1d{-300}; + const Int_t eb_TrigWinMaxTrd1d{300}; + const Int_t eb_TrigWinMinTrd2d{-200}; + const Int_t eb_TrigWinMaxTrd2d{200}; + const Int_t eb_TrigWinMinTof{-20}; + const Int_t eb_TrigWinMaxTof{60}; + const Int_t eb_TrigWinMinRich{-60}; + const Int_t eb_TrigWinMaxRich{60}; + + if (bEvB) { + CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents(); + evBuildRaw->SetVerbose(3); + //Choose between NoOverlap, MergeOverlap, AllowOverlap + evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap); + //For time being it is needed. We will remove CbmMuchBeamTimeDigi. + evBuildRaw->ChangeMuchBeamtimeDigiFlag(); + // Remove detectors where digis not found + if (!bRICH || !geoSetup->IsActive(ECbmModuleId::kRich)) evBuildRaw->RemoveDetector(kRawEventBuilderDetRich); + if (!bMUCH || !geoSetup->IsActive(ECbmModuleId::kMuch)) evBuildRaw->RemoveDetector(kRawEventBuilderDetMuch); + if (!bPSD || !geoSetup->IsActive(ECbmModuleId::kPsd)) evBuildRaw->RemoveDetector(kRawEventBuilderDetPsd); + if (!bTRD || !geoSetup->IsActive(ECbmModuleId::kTrd)) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd); + //TODO temporary exclude 2D from event building (asked by M. Beyer) + if (!bTRD2d) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd2D); + if (!bSTS || !geoSetup->IsActive(ECbmModuleId::kSts)) evBuildRaw->RemoveDetector(kRawEventBuilderDetSts); + if (!bTOF || !geoSetup->IsActive(ECbmModuleId::kTof)) evBuildRaw->RemoveDetector(kRawEventBuilderDetTof); + if (!bFSD || !geoSetup->IsActive(ECbmModuleId::kFsd)) evBuildRaw->RemoveDetector(kRawEventBuilderDetFsd); + + // Set Bmon as reference detector. + // Select only Bmon1 [newly installed for mcbm2024 data] (AB) + evBuildRaw->SetReferenceDetector(kRawEventBuilderDetBmon, {0, 1}); + // For making MuCh as seed detector + // evBuildRaw->SetReferenceDetector(kRawEventBuilderDetMuch); + + // Use sliding window seed builder with STS + //evBuildRaw->SetReferenceDetector(kRawEventBuilderDetUndef); + //evBuildRaw->AddSeedTimeFillerToList(kRawEventBuilderDetTof); + //evBuildRaw->SetSlidingWindowSeedFinder(10, 10, 50); + //evBuildRaw->SetSeedFinderQa(true); // optional QA information for seed finder + //evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, 1000); + //evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1); + //evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, -500, 500); + + // void SetTsParameters(double TsStartTime, double TsLength, double TsOverLength): + // => TsStartTime=0, TsLength=128 + 1.28 ms in 2022, TsOverLength=1.28 ms (1MS) in mCBM2022 + evBuildRaw->SetTsParameters(0.0, 1.28e8, 1.28e6); + + if (geoSetup->IsActive(ECbmModuleId::kTof)) { + evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kTof, eb_TriggerMinNumberTof); + evBuildRaw->SetTriggerMinLayersNumber(ECbmModuleId::kTof, eb_TriggerMinNumberTofLayers); + evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kTof, -1); + } + + evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kBmon, eb_TriggerMinNumberBmon); + evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kBmon, eb_TriggerMaxNumberBmon); + + if (geoSetup->IsActive(ECbmModuleId::kSts)) { + evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, eb_TriggerMinNumberSts); + evBuildRaw->SetTriggerMinLayersNumber(ECbmModuleId::kSts, eb_TriggerMinNumberStsLayers); + evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1); + } + if (bMUCH && geoSetup->IsActive(ECbmModuleId::kMuch)) { + evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kMuch, eb_TriggerMinNumberMuch); + evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kMuch, -1); + } + if (geoSetup->IsActive(ECbmModuleId::kRich)) { + evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kRich, eb_TriggerMinNumberRich); + evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kRich, -1); + } + + evBuildRaw->SetTriggerWindow(ECbmModuleId::kBmon, eb_TrigWinMinBmon, eb_TrigWinMaxBmon); + if (geoSetup->IsActive(ECbmModuleId::kTof)) + evBuildRaw->SetTriggerWindow(ECbmModuleId::kTof, eb_TrigWinMinTof, eb_TrigWinMaxTof); + if (geoSetup->IsActive(ECbmModuleId::kSts)) + evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, eb_TrigWinMinSts, eb_TrigWinMaxSts); + if (bMUCH && geoSetup->IsActive(ECbmModuleId::kMuch)) + evBuildRaw->SetTriggerWindow(ECbmModuleId::kMuch, eb_TrigWinMinMuch, eb_TrigWinMaxMuch); + if (geoSetup->IsActive(ECbmModuleId::kTrd)) + evBuildRaw->SetTriggerWindow(ECbmModuleId::kTrd, eb_TrigWinMinTrd1d, eb_TrigWinMaxTrd1d); + if (geoSetup->IsActive(ECbmModuleId::kTrd)) + evBuildRaw->SetTriggerWindow(ECbmModuleId::kTrd2d, eb_TrigWinMinTrd2d, eb_TrigWinMaxTrd2d); + if (geoSetup->IsActive(ECbmModuleId::kRich)) + evBuildRaw->SetTriggerWindow(ECbmModuleId::kRich, eb_TrigWinMinRich, eb_TrigWinMaxRich); + + run->AddTask(evBuildRaw); + } + // ------------------------------------------------------------------------ + + + // ----- Reconstruction tasks ----------------------------------------- + + + // ========================================================================= + // === local STS Reconstruction === + // ========================================================================= + + if (bSTS && geoSetup->IsActive(ECbmModuleId::kSts)) { + CbmRecoSts* recoSts = new CbmRecoSts(); + if (bEvB) { + recoSts->SetMode(ECbmRecoMode::EventByEvent); + } + else { + recoSts->SetMode(ECbmRecoMode::Timeslice); + } + recoSts->SetVerbose(3); + recoSts->SetTimeCutDigisAbs(20.0); // cluster finder: time cut in ns + recoSts->SetTimeCutClustersAbs(20.0); // hit finder: time cut in ns + + // Sensor params + CbmStsParSensor sensor6cm(CbmStsSensorClass::kDssdStereo); + sensor6cm.SetPar(0, 6.2092); // Extension in x + sensor6cm.SetPar(1, 6.2); // Extension in y + sensor6cm.SetPar(2, 0.03); // Extension in z + sensor6cm.SetPar(3, 5.9692); // Active size in y + sensor6cm.SetPar(4, 1024.); // Number of strips front side + sensor6cm.SetPar(5, 1024.); // Number of strips back side + sensor6cm.SetPar(6, 0.0058); // Strip pitch front side + sensor6cm.SetPar(7, 0.0058); // Strip pitch back side + sensor6cm.SetPar(8, 0.0); // Stereo angle front side + sensor6cm.SetPar(9, 7.5); // Stereo angle back side + + CbmStsParSensor sensor12cm(sensor6cm); // copy all parameters, change then only the y size + sensor12cm.SetPar(1, 12.4); // Extension in y + sensor12cm.SetPar(3, 12.1692); // Active size in y + + // --- Addresses for sensors + // --- They are defined in each station as sensor 1, module 1, halfladderD (2), ladder 1 + // Int_t GetAddress(UInt_t unit = 0, UInt_t ladder = 0, UInt_t halfladder = 0, UInt_t module = 0, UInt_t sensor = 0, + // UInt_t side = 0, UInt_t version = kCurrentVersion); + + // setup available starting with sts_v24b. Available if reconstructing from rra files (AB 14.05.2024) + Int_t stsAddress00 = CbmStsAddress::GetAddress(0, 0, 0, 0, 0, 0); // U0 L0 M0 6 cm + Int_t stsAddress01 = CbmStsAddress::GetAddress(1, 0, 1, 0, 0, 0); // U1 L0 M0 6 cm + Int_t stsAddress02 = CbmStsAddress::GetAddress(1, 0, 1, 1, 0, 0); // U1 L0 M1 6 cm + Int_t stsAddress03 = CbmStsAddress::GetAddress(1, 1, 1, 0, 0, 0); // U1 L1 M0 6 cm + Int_t stsAddress04 = CbmStsAddress::GetAddress(1, 1, 1, 1, 0, 0); // U1 L1 M1 6 cm + Int_t stsAddress05 = CbmStsAddress::GetAddress(2, 0, 1, 0, 0, 0); // U2 L0 M0 6 cm + Int_t stsAddress06 = CbmStsAddress::GetAddress(2, 0, 1, 1, 0, 0); // U2 L0 M1 12 cm + Int_t stsAddress07 = CbmStsAddress::GetAddress(2, 1, 1, 0, 0, 0); // U2 L1 M0 6 cm + Int_t stsAddress08 = CbmStsAddress::GetAddress(2, 1, 1, 1, 0, 0); // U2 L1 M1 12 cm + Int_t stsAddress09 = CbmStsAddress::GetAddress(2, 2, 1, 0, 0, 0); // U2 L2 M0 6 cm + Int_t stsAddress10 = CbmStsAddress::GetAddress(2, 2, 1, 1, 0, 0); // U2 L2 M1 6 cm + Int_t stsAddress11 = CbmStsAddress::GetAddress(2, 2, 1, 2, 0, 0); // U2 L2 M2 6 cm + + std::cout << "STS address00 " << std::dec << stsAddress00 << " " << std::hex << stsAddress00 << std::endl; + std::cout << "STS address01 " << std::dec << stsAddress01 << " " << std::hex << stsAddress01 << std::endl; + std::cout << "STS address02 " << std::dec << stsAddress02 << " " << std::hex << stsAddress02 << std::endl; + std::cout << "STS address03 " << std::dec << stsAddress03 << " " << std::hex << stsAddress03 << std::endl; + std::cout << "STS address04 " << std::dec << stsAddress04 << " " << std::hex << stsAddress04 << std::endl; + std::cout << "STS address05 " << std::dec << stsAddress05 << " " << std::hex << stsAddress05 << std::endl; + std::cout << "STS address06 " << std::dec << stsAddress06 << " " << std::hex << stsAddress06 << std::endl; + std::cout << "STS address07 " << std::dec << stsAddress07 << " " << std::hex << stsAddress07 << std::endl; + std::cout << "STS address08 " << std::dec << stsAddress08 << " " << std::hex << stsAddress08 << std::endl; + std::cout << "STS address09 " << std::dec << stsAddress09 << " " << std::hex << stsAddress09 << std::endl; + std::cout << "STS address10 " << std::dec << stsAddress10 << " " << std::hex << stsAddress10 << std::endl; + std::cout << "STS address11 " << std::dec << stsAddress11 << " " << std::hex << stsAddress11 << std::endl; + + // --- Now we can define the sensor parameter set and tell recoSts to use it + auto sensorParSet = new CbmStsParSetSensor("CbmStsParSetSensor", "STS sensor parameters" + "mcbm2024"); + if (stsAddress00) sensorParSet->SetParSensor(stsAddress00, sensor6cm); + sensorParSet->SetParSensor(stsAddress01, sensor6cm); + sensorParSet->SetParSensor(stsAddress02, sensor6cm); + sensorParSet->SetParSensor(stsAddress03, sensor6cm); + sensorParSet->SetParSensor(stsAddress04, sensor6cm); + sensorParSet->SetParSensor(stsAddress05, sensor6cm); + sensorParSet->SetParSensor(stsAddress06, sensor12cm); + sensorParSet->SetParSensor(stsAddress07, sensor6cm); + sensorParSet->SetParSensor(stsAddress08, sensor12cm); + sensorParSet->SetParSensor(stsAddress09, sensor6cm); + sensorParSet->SetParSensor(stsAddress10, sensor6cm); + sensorParSet->SetParSensor(stsAddress11, sensor6cm); + recoSts->UseSensorParSet(sensorParSet); + + // ASIC params: #ADC channels, dyn. range, threshold, time resol., dead time, + // noise RMS, zero-threshold crossing rate + auto parAsic = new CbmStsParAsic(128, 31, 75000., 3000., 5., 800., 1000., 3.9789e-3); + + // Module params: number of channels, number of channels per ASIC + auto parMod = new CbmStsParModule(2048, 128); + parMod->SetAllAsics(*parAsic); + recoSts->UseModulePar(parMod); + + // Sensor conditions: full depletion voltage, bias voltage, temperature, + // coupling capacitance, inter-strip capacitance + auto sensorCond = new CbmStsParSensorCond(70., 140., 268., 17.5, 1.); + recoSts->UseSensorCond(sensorCond); + + run->AddTask(recoSts); + std::cout << "-I- : Added task " << recoSts->GetName() << std::endl; + // ------------------------------------------------------------------------ + } + + // ========================================================================= + // === local MUCH Reconstruction === + // ========================================================================= + if (bMUCH) { + TString muchGeoTag; + if (geoSetup->GetGeoTag(ECbmModuleId::kMuch, muchGeoTag)) { + // --- Parameter file name + TString geoTag; + geoSetup->GetGeoTag(ECbmModuleId::kMuch, geoTag); + Int_t muchFlag = 0; + if (geoTag.Contains("mcbm")) muchFlag = 1; + TString sectorFile = gSystem->Getenv("VMCWORKDIR"); + sectorFile += "/parameters/much/much_" + geoTag(0, 4) + "_mcbm_digi_sector.root"; + //sectorFile += "/parameters/much/much_v22j_mcbm_digi_sector.root"; + std::cout << "Using parameter file " << sectorFile << std::endl; + + // --- Initialization of the digi scheme + auto muchGeoScheme = CbmMuchGeoScheme::Instance(); + if (!muchGeoScheme->IsInitialized()) { + muchGeoScheme->Init(sectorFile.Data(), muchFlag); + } + // --- Hit finder for GEMs + FairTask* muchReco = new CbmMuchFindHitsGem(sectorFile.Data(), muchFlag); + run->AddTask(muchReco); + std::cout << "-I- " << myName << ": Added task " << muchReco->GetName() << " with parameter file " << sectorFile + << std::endl; + } + } + // ------------------------------------------------------------------------ + + // ========================================================================= + // === local TRD Reconstruction === + // ========================================================================= + + if (bTRD && geoSetup->IsActive(ECbmModuleId::kTrd)) { + CbmTrdClusterFinder* trdCluster; + Double_t triggerThreshold = 0.5e-6; // SIS100 + + trdCluster = new CbmTrdClusterFinder(); + trdCluster->SetNeighbourEnable(true, false); + trdCluster->SetMinimumChargeTH(triggerThreshold); + trdCluster->SetRowMerger(true); + CbmTrdClusterFinder::UseRecoHelpers(); + run->AddTask(trdCluster); + std::cout << "-I- : Added task " << trdCluster->GetName() << std::endl; + + CbmTrdHitProducer* trdHit = new CbmTrdHitProducer(); + //trdHit->SetHitTimeOffset(363); // hit time synchronization for TRD2D determined on run 2391 + trdHit->SetHitTimeOffset(80); // hit time synchronization for TRD2D determined on run 2914 + run->AddTask(trdHit); + std::cout << "-I- : Added task " << trdHit->GetName() << std::endl; + } + + // ========================================================================= + // === RICH Reconstruction === + // ========================================================================= + + if (bRICH && geoSetup->IsActive(ECbmModuleId::kRich)) { + // ----- Local reconstruction of RICH Hits ------------------------------ + CbmRichMCbmHitProducer* hitProd = new CbmRichMCbmHitProducer(); + hitProd->SetMappingFile(std::string(srcDir.Data()) + + "/macro/rich/mcbm/beamtime/mRICH_Mapping_vert_20190318_elView.geo"); + hitProd->SetIcdFilenameBase(std::string(srcDir.Data()) + "/macro/beamtime/mcbm2022/icd_offset_it"); + hitProd->setToTLimits(23.7, 30.0); + hitProd->applyToTCut(); + hitProd->applyICDCorrection(); + run->AddTask(hitProd); + // ------------------------------------------------------------------------ + + // ----- Local reconstruction in RICh -> Finding of Rings --------------- + CbmRichReconstruction* richReco = new CbmRichReconstruction(); + richReco->UseMCbmSetup(); + run->AddTask(richReco); + // ------------------------------------------------------------------------ + } + + // ========================================================================= + // === TOF Hitfinding === + // ========================================================================= + + TString parPath = srcDir + "/parameters/mcbm/"; + if (bTOF && geoSetup->IsActive(ECbmModuleId::kTof)) { + TString cFname; + switch (iTofCluMode) { + case 1: { + // ----- TOF defaults ------------------------ + Int_t calMode = 93; + Int_t calSel = 1; + Int_t calSm = 0; + Int_t RefSel = 0; + Double_t dDeadtime = 50.; + Int_t iSel2 = 500; + Bool_t bOut = kTRUE; + + // ------------------------------------------------------------------------ + gROOT->LoadMacro(srcDir + "/macro/beamtime/mcbm2022/ini_tof_clusterizer.C"); + Char_t* cCmd = + Form("ini_tof_clusterizer(%d,%d,%d,%d,\"%s\",%d,%d,%d,%f,\"%s\",\"%s\")", calMode, calSel, calSm, RefSel, + cFileId.Data(), iCalSet, (Int_t) bOut, iSel2, dDeadtime, cCalId.Data(), parPath.Data()); + + cout << "<I> " << cCmd << endl; + gInterpreter->ProcessLine(cCmd); + // disable histogramming + CbmTofEventClusterizer* tofClust = CbmTofEventClusterizer::Instance(); + tofClust->SetDutId(-1); // to disable histogramming + } break; + + default: { + ; + } + } + // ------------------------------------------------------------------------- + + // ========================================================================= + // === Tof Tracking === + // ========================================================================= + + if (bTOFtr) { + cout << "<I> Initialize Tof tracker by ini_trks" << endl; + + gROOT->LoadMacro(srcDir + "/macro/beamtime/mcbm2022/ini_tof_trks.C"); + Char_t* cCmd = + Form("ini_tof_trks(%d,%d,%d,%6.2f,%8.1f,\"%s\",%d,%d,%d,%f,\"%s\")", iSel, iTrackingSetup, iGenCor, dScalFac, + dChi2Lim2, cCalId.Data(), (Int_t) bUseSigCalib, iCalOpt, iTrkPar, dTOffScal, parPath.Data()); + cout << "<I> " << cCmd << endl; + gInterpreter->ProcessLine(cCmd); + + CbmTofFindTracks* tofFindTracks = CbmTofFindTracks::Instance(); + Int_t iNStations = tofFindTracks->GetNStations(); + } + } + + + //Constant Field + // CbmConstField *fMagField = new CbmConstField(); + // fMagField->SetFieldXYZ(0, 0 ,0 ); // values are in kG + // fMagField->SetFieldRegions(-10, -10 ,-10 , 10, 10 , 10 ); + // run->SetField(fMagField); + + // ========================================================================= + // === L1 === + // ========================================================================= + if (bL1) { + run->AddTask(new CbmTrackingDetectorInterfaceInit()); + + CbmL1* l1 = new CbmL1("L1"); + l1->SetMcbmMode(); + l1->SetVerbose(3); + run->AddTask(l1); + + CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder(); + FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder); + run->AddTask(globalFindTracks); + + if (bPV) { + CbmKF* KF = new CbmKF(); + run->AddTask(KF); + + auto* pvFinder = new CbmPVFinderKFGlobal(); + CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder); + findVertex->SetTrackType(ECbmDataType::kGlobalTrack); + run->AddTask(findVertex); + } + } + // ========================================================================= + // === QA === + // ========================================================================= + if (bQA) { + CbmRecoQaTask* recoQa = new CbmRecoQaTask(); + recoQa->SetSetupClass(CbmRecoQaTask::kMcbm24); + // USER : Uncomment this line if you like to select only track multiplicities 1 and 2 per event + // recoQa->AddEventFilter(CbmRecoQaTask::EventFilter::eEventCut::kMultTrk)->SetFilter({1, 2}); + // USER : All track cuts are EXCLUSIVE (they ALL have to be met in order that tracks are selected) + // USER : Uncomment this line if you like to select ONLY tracks with AT LEAST 3 STS hits + // recoQa->AddTrackFilter(CbmRecoQaTask::TrackFilter::eTrackCut::kSts)->SetFilter({3}); + // USER : Uncomment this line if you like to select ONLY tracks with AT LEAST 2 TRD hits + // recoQa->AddTrackFilter(CbmRecoQaTask::TrackFilter::eTrackCut::kTrd)->SetFilter({2}); + // USER : Uncomment this line if you like to select ONLY tracks with AT LEAST 1 ToF hit + // recoQa->AddTrackFilter(CbmRecoQaTask::TrackFilter::eTrackCut::kTof)->SetFilter({1}); + run->AddTask(recoQa); + } + // ------------------------------------------------------------------------ + + // ----- 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(); + FairParRootFileIo* parIo3 = new FairParRootFileIo(); + //parIo1->open(parFileIn.Data(), "READ"); + //rtdb->setFirstInput(parIo1); + parIo2->open(parFileList, "in"); + rtdb->setSecondInput(parIo2); + parIo3->open(parFileOut.Data(), "RECREATE"); + // ------------------------------------------------------------------------ + rtdb->setOutput(parIo3); + + // ----- Run initialisation ------------------------------------------- + std::cout << std::endl; + std::cout << "-I- " << myName << ": Initialise run" << std::endl; + run->Init(); + rtdb->print(); + // ------------------------------------------------------------------------ + + + // ----- Start run ---------------------------------------------------- + std::cout << std::endl << std::endl; + std::cout << "-I- " << myName << ": Starting run" << std::endl; + run->Run(0, nTimeslices); + // ------------------------------------------------------------------------ + + + // ----- Finish ------------------------------------------------------- + + rtdb->print(); + rtdb->saveOutput(); + run->CreateGeometryFile(geoFileOut); + + 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 " << parFileOut << std::endl; + std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << std::endl; + std::cout << std::endl; + // ------------------------------------------------------------------------ + + + // ----- Resource monitoring ------------------------------------------ + // Extract the maximal used memory an add is as Dart measurement + // This line is filtered by CTest and the value send to CDash + 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; + // ------------------------------------------------------------------------ + + + // ----- Function needed for CTest runtime dependency ----------------- + RemoveGeoManager(); + // ------------------------------------------------------------------------ + + /// --- Screen output for automatic tests + std::cout << " Test passed" << std::endl; + std::cout << " All ok " << std::endl; + + return kTRUE; +}