Skip to content
Snippets Groups Projects
Commit d75b7b39 authored by Sergei Zharko's avatar Sergei Zharko
Browse files

mcbm_event_reco_L1.C for mcbm2025

parent fbaf37cb
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
macro/beamtime/common/CMakeLists.txt
+ 9
− 0
View file @ d75b7b39
...@@ -24,3 +24,12 @@ Install(DIRECTORY online ...@@ -24,3 +24,12 @@ Install(DIRECTORY online
Install(DIRECTORY data Install(DIRECTORY data
DESTINATION share/cbmroot/macro/beamtime/common DESTINATION share/cbmroot/macro/beamtime/common
PATTERN "*" EXCLUDE) 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} )
macro/beamtime/mcbm2025/mcbm_event_reco_L1.C 0 → 100644
+ 738
− 0
View file @ d75b7b39
/* 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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment