From d547d69a5142a1da1dd3c108c064498f879ace5e Mon Sep 17 00:00:00 2001
From: Adrian Weber <adrian.a.weber@physik.uni-giessen.de>
Date: Wed, 16 Feb 2022 11:14:30 +0100
Subject: [PATCH] add data reconstruction macro for 2021, analysis macro for
alignment + alignment matrix generator
---
.../beamtime/mcbm2021/create_alignment_file.C | 67 +
macro/beamtime/mcbm2021/mcbm_event_reco.C | 969 +++++++++++++
.../mcbm2021/mcbm_inspect_alignment.C | 1214 +++++++++++++++++
reco/detectors/tof/CbmTofFindTracks.cxx | 1 +
4 files changed, 2251 insertions(+)
create mode 100644 macro/beamtime/mcbm2021/create_alignment_file.C
create mode 100644 macro/beamtime/mcbm2021/mcbm_event_reco.C
create mode 100644 macro/beamtime/mcbm2021/mcbm_inspect_alignment.C
diff --git a/macro/beamtime/mcbm2021/create_alignment_file.C b/macro/beamtime/mcbm2021/create_alignment_file.C
new file mode 100644
index 0000000000..72803fcacc
--- /dev/null
+++ b/macro/beamtime/mcbm2021/create_alignment_file.C
@@ -0,0 +1,67 @@
+/* Copyright (C) 2022 UGiessen, Giessen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Florian Uhlig, Adrian Weber [committer]*/
+
+#include <TFile.h>
+#include <TGeoMatrix.h>
+
+#include <map>
+#include <string>
+
+std::pair<std::string, TGeoHMatrix> AlignNode(std::string path, double shiftX, double shiftY, double shiftZ,
+ double rotX, double rotY, double rotZ)
+{
+
+ TGeoHMatrix result;
+ result.SetDx(shiftX);
+ result.SetDy(shiftY);
+ result.SetDz(shiftZ);
+ result.RotateX(rotX);
+ result.RotateY(rotY);
+ result.RotateZ(rotZ);
+
+ std::cout << "Alignment matrix for node " << path << " is: " << std::endl;
+ result.Print();
+ std::cout << std::endl;
+
+ return std::pair<std::string, TGeoHMatrix>(path, result);
+}
+
+
+int create_alignment_file()
+{
+ // 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;
+
+ // ---------------- STS ----------------------------//
+ // Align full STS
+ matrices.insert(AlignNode("/cave_1/sts_v22c_mcbm_0", 0.275, -0.89, -.5, 0., 0., 0.));
+
+ // Align individual STS Units
+ // Station 1
+ // Unit 0
+ matrices.insert(AlignNode("/cave_1/sts_v22c_mcbm_0/Station01_1/Ladder09_1", 0.0, 0.045, 0., 0., 0., 0.));
+ // Unit 1
+ matrices.insert(AlignNode("/cave_1/sts_v22c_mcbm_0/Station01_1/Ladder09_2", -0.04, 0.06, 0., 0., 0., 0.));
+
+ // Station 2
+ // Unit 2
+ matrices.insert(AlignNode("/cave_1/sts_v22c_mcbm_0/Station02_2/Ladder10_2", 0.0, -0.11, 0., 0., 0., 0.));
+
+ // ---------------- TOF ----------------------------//
+ // Align full Tof
+ matrices.insert(AlignNode("/cave_1/tof_v21d_mcbm_0/tof_v21d_mcbmStand_1", 0.85, -1.05, 0.0, 0., 0.5, 0.));
+
+
+ // save matrices to disk
+ TFile* misalignmentMatrixRootfile = new TFile("AlignmentMatrices.root", "RECREATE");
+ if (misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->WriteObject(&matrices, "MisalignMatrices");
+ misalignmentMatrixRootfile->Write();
+ misalignmentMatrixRootfile->Close();
+ }
+
+ return 0;
+}
diff --git a/macro/beamtime/mcbm2021/mcbm_event_reco.C b/macro/beamtime/mcbm2021/mcbm_event_reco.C
new file mode 100644
index 0000000000..d44e297085
--- /dev/null
+++ b/macro/beamtime/mcbm2021/mcbm_event_reco.C
@@ -0,0 +1,969 @@
+/* Copyright (C) 2020 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau, Adrian Weber [committer] */
+
+// --------------------------------------------------------------------------
+//
+// Macro for reconstruction of mcbm data (2021)
+// Combined reconstruction (cluster + hit finder) for different subsystems.
+//
+// --------------------------------------------------------------------------
+
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+/// FIXME: Disable clang formatting to keep easy parameters overview
+/* clang-format off */
+Bool_t mcbm_event_reco(UInt_t uRunId = 1588,
+ Int_t nTimeslices = 20,
+ TString sInpDir = "/data/cbmroot/files/mTofCriPar2",
+ TString sOutDir = "rec/tofPar2_noAlign_L1/",
+ TString alignmentMatrixFileName = "AlignmentMatrices.root",
+ Int_t iUnpFileIndex = -1)
+{
+ /// FIXME: Re-enable clang formatting after parameters initial values setting
+ /* clang-format on */
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "INFO";
+ TString logVerbosity = "LOW";
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ TString myName = "mcbm_event_reco"; // 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);
+
+ /// Initial pattern
+ TString inFile = sInpDir + "/" + sRunId + ".digi";
+
+ //TString parFileIn = sInpDir + "/unp_mcbm_params_" + sRunId;
+ TString parFileOut = sOutDir + "/reco_event_mcbm_params_" + sRunId;
+ TString outFile = sOutDir + "/reco_event_mcbm_" + sRunId;
+
+ // Your folder with the Tof Calibration files;
+ TString TofFileFolder = "/data/cbmroot/files/tofCal/mTofCriPar2/";
+
+
+ /// 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!
+ parFileOut += TString::Format("_%02u", iUnpFileIndex);
+ outFile += TString::Format("_%02u", iUnpFileIndex);
+ } // if ( 0 <= uUnpFileIndex )
+ /// Add ROOT file suffix
+ inFile += ".root";
+ // parFileIn += ".root";
+ parFileOut += ".root";
+ outFile += ".root";
+ // ---------------------------------------------
+
+ // ----- EventBuilder Settings----------------
+ const Int_t eb_TriggerMinNumberT0 {0};
+ const Int_t eb_TriggerMinNumberSts {0};
+ const Int_t eb_TriggerMinNumberMuch {0};
+ const Int_t eb_TriggerMinNumberTof {2};
+ const Int_t eb_TriggerMinNumberRich {5};
+
+ // ----- 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 = 20; //500;
+
+ // Tracking
+ Int_t iSel = 1; //500;//910041;
+ Int_t iTrackingSetup = 2;
+ Int_t iGenCor = 1;
+ Double_t dScalFac = 1.;
+ Double_t dChi2Lim2 = 500.;
+ Bool_t bUseSigCalib = kFALSE;
+ Int_t iCalOpt = 0;
+ Int_t iTrkPar = 3;
+ // ------------------------------------------------------------------------
+
+ // ----- 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";
+ Int_t iCalSet = 30040500; // calibration settings
+ if (uRunId >= 759) iCalSet = 10020500;
+ if (uRunId >= 812) iCalSet = 10020500;
+ if (uRunId >= 1588) iCalSet = 12002002;
+
+ 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)
+ TString geoSetupTag = "mcbm_beam_2021_07_surveyed";
+ TString geoFile = srcDir + "/macro/mcbm/data/" + geoSetupTag + ".geo.root";
+ CbmSetup* geoSetup = CbmSetup::Instance();
+ geoSetup->LoadSetup(geoSetupTag);
+
+ // You can modify the pre-defined setup by using
+ geoSetup->SetActive(ECbmModuleId::kMvd, kFALSE);
+ geoSetup->SetActive(ECbmModuleId::kSts, kTRUE);
+ geoSetup->SetActive(ECbmModuleId::kMuch, kFALSE);
+ geoSetup->SetActive(ECbmModuleId::kRich, kFALSE);
+ geoSetup->SetActive(ECbmModuleId::kTrd, kTRUE);
+ geoSetup->SetActive(ECbmModuleId::kTof, kTRUE);
+ geoSetup->SetActive(ECbmModuleId::kPsd, kFALSE);
+
+ //----- 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 geoTag;
+ if (geoSetup->IsActive(ECbmModuleId::kTof)) {
+ geoSetup->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;
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- 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("reco", "reco.monitor");
+ 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 (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---------------------------------------
+ CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents();
+
+ //Choose between NoOverlap, MergeOverlap, AllowOverlap
+ evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap);
+
+ // Remove detectors where digis not found
+ if (!geoSetup->IsActive(ECbmModuleId::kRich)) evBuildRaw->RemoveDetector(kRawEventBuilderDetRich);
+ if (!geoSetup->IsActive(ECbmModuleId::kMuch)) evBuildRaw->RemoveDetector(kRawEventBuilderDetMuch);
+ if (!geoSetup->IsActive(ECbmModuleId::kPsd)) evBuildRaw->RemoveDetector(kRawEventBuilderDetPsd);
+ if (!geoSetup->IsActive(ECbmModuleId::kTrd)) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd);
+ if (!geoSetup->IsActive(ECbmModuleId::kSts)) evBuildRaw->RemoveDetector(kRawEventBuilderDetSts);
+ if (!geoSetup->IsActive(ECbmModuleId::kTof)) evBuildRaw->RemoveDetector(kRawEventBuilderDetTof);
+
+ // Set TOF as reference detector
+ evBuildRaw->SetReferenceDetector(kRawEventBuilderDetTof);
+
+ // void SetTsParameters(double TsStartTime, double TsLength, double TsOverLength): TsStartTime=0, TsLength=256ms in 2021, TsOverLength=TS overlap, not used in mCBM2021
+ evBuildRaw->SetTsParameters(0.0, 2.56e8, 0.0);
+
+ if (geoSetup->IsActive(ECbmModuleId::kTof))
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kTof, eb_TriggerMinNumberTof);
+
+ if (geoSetup->IsActive(ECbmModuleId::kTof)) evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kTof, -1);
+
+ //evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, eb_TriggerMinNumberSts);
+ //evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1);
+
+ //evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kRich, eb_TriggerMinNumberRich);
+ //evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kRich, -1);
+
+ if (geoSetup->IsActive(ECbmModuleId::kTof)) evBuildRaw->SetTriggerWindow(ECbmModuleId::kTof, -50, 50);
+ if (geoSetup->IsActive(ECbmModuleId::kSts)) evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, -50, 50);
+ if (geoSetup->IsActive(ECbmModuleId::kTrd)) evBuildRaw->SetTriggerWindow(ECbmModuleId::kTrd, -200, 200);
+ if (geoSetup->IsActive(ECbmModuleId::kRich)) evBuildRaw->SetTriggerWindow(ECbmModuleId::kRich, -50, 50);
+
+ run->AddTask(evBuildRaw);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Reconstruction tasks -----------------------------------------
+
+
+ // =========================================================================
+ // === local STS Reconstruction ===
+ // =========================================================================
+
+ if (geoSetup->IsActive(ECbmModuleId::kSts)) {
+ CbmRecoSts* recoSts = new CbmRecoSts();
+ recoSts->SetMode(kCbmRecoEvent);
+
+ 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);
+
+ Int_t stsAddress01 = CbmStsAddress::GetAddress(0, 0, 1, 0, 0, 0); // U0 L0 M0 6 cm
+ Int_t stsAddress02 = CbmStsAddress::GetAddress(0, 0, 1, 1, 0, 0); // U0 L0 M1 6 cm
+ Int_t stsAddress03 = CbmStsAddress::GetAddress(0, 1, 1, 0, 0, 0); // U0 L1 M0 6 cm
+ Int_t stsAddress04 = CbmStsAddress::GetAddress(0, 1, 1, 1, 0, 0); // U0 L1 M1 6 cm
+ Int_t stsAddress05 = CbmStsAddress::GetAddress(1, 0, 1, 0, 0, 0); // U1 L0 M0 6 cm
+ Int_t stsAddress06 = CbmStsAddress::GetAddress(1, 0, 1, 1, 0, 0); // U1 L0 M1 12 cm
+ Int_t stsAddress07 = CbmStsAddress::GetAddress(1, 1, 1, 0, 0, 0); // U1 L1 M0 6 cm
+ Int_t stsAddress08 = CbmStsAddress::GetAddress(1, 1, 1, 1, 0, 0); // U1 L1 M1 12 cm
+ Int_t stsAddress09 = CbmStsAddress::GetAddress(1, 2, 1, 0, 0, 0); // U1 L2 M0 6 cm
+ Int_t stsAddress10 = CbmStsAddress::GetAddress(1, 2, 1, 1, 0, 0); // U1 L2 M1 6 cm
+ Int_t stsAddress11 = CbmStsAddress::GetAddress(1, 2, 1, 2, 0, 0); // U1 L2 M2 6 cm
+
+
+ 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"
+ "mcbm2021");
+ 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, 32, 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 TRD Reconstruction ===
+ // =========================================================================
+
+ CbmTrdClusterFinder* trdCluster;
+ if (geoSetup->IsActive(ECbmModuleId::kTrd)) {
+ Double_t triggerThreshold = 0.5e-6; // SIS100
+
+ trdCluster = new CbmTrdClusterFinder();
+ trdCluster->SetNeighbourEnable(true, false);
+ trdCluster->SetMinimumChargeTH(triggerThreshold);
+ trdCluster->SetRowMerger(true);
+ run->AddTask(trdCluster);
+ std::cout << "-I- : Added task " << trdCluster->GetName() << std::endl;
+
+ CbmTrdHitProducer* trdHit = new CbmTrdHitProducer();
+ run->AddTask(trdHit);
+ std::cout << "-I- : Added task " << trdHit->GetName() << std::endl;
+ }
+
+
+ // =========================================================================
+ // === RICH Reconstruction ===
+ // =========================================================================
+
+ if (geoSetup->IsActive(ECbmModuleId::kRich)) {
+ // ----- Local reconstruction of RICH Hits ------------------------------
+ CbmRichMCbmHitProducer* hitProd = new CbmRichMCbmHitProducer();
+ hitProd->setToTLimits(23.7, 30.0);
+ hitProd->applyToTCut();
+ run->AddTask(hitProd);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Local reconstruction in RICh -> Finding of Rings ---------------
+ CbmRichReconstruction* richReco = new CbmRichReconstruction();
+ richReco->UseMCbmSetup();
+ run->AddTask(richReco);
+ // ------------------------------------------------------------------------
+ }
+
+ // =========================================================================
+ // === TOF Hitfinding ===
+ // =========================================================================
+
+ if (geoSetup->IsActive(ECbmModuleId::kTof)) {
+ TString cFname;
+ switch (iTofCluMode) {
+ case 1: {
+ CbmTofEventClusterizer* tofCluster = new CbmTofEventClusterizer("TOF Event Clusterizer", 0, 1);
+ cFname = Form("/%s_set%09d_%02d_%01dtofClust.hst.root", cCalId.Data(), iCalSet, calMode, calSel);
+ tofCluster->SetCalParFileName(TofFileFolder + cFname);
+ tofCluster->SetCalMode(calMode);
+ tofCluster->SetCalSel(calSel);
+ tofCluster->SetCaldXdYMax(300.); // geometrical matching window in cm
+ tofCluster->SetCalCluMulMax(3.); // Max Counter Cluster Multiplicity for filling calib histos
+ tofCluster->SetCalRpc(calSm); // select detector for calibration update
+ tofCluster->SetTRefId(RefSel); // reference trigger for offset calculation
+ tofCluster->SetTotMax(20.); // Tot upper limit for walk corection
+ tofCluster->SetTotMin(0.); //(12000.); // Tot lower limit for walk correction
+ tofCluster->SetTotPreRange(5.); // effective lower Tot limit in ns from peak position
+ tofCluster->SetTotMean(5.); // Tot calibration target value in ns
+ tofCluster->SetMaxTimeDist(1.0); // default cluster range in ns
+ tofCluster->SetDelTofMax(50.); // acceptance range for cluster distance in ns (!)
+ tofCluster->SetSel2MulMax(3); // limit Multiplicity in 2nd selector
+ tofCluster->SetChannelDeadtime(dDeadtime); // artificial deadtime in ns
+ tofCluster->SetEnableAvWalk(kFALSE);
+ //tofCluster->SetEnableMatchPosScaling(kFALSE); // turn off projection to nominal target
+ tofCluster->SetYFitMin(1.E4);
+ tofCluster->SetToDAv(0.04);
+ tofCluster->SetIdMode(1); // calibrate on module level
+ tofCluster->SetTRefDifMax(2.0); // in ns
+ tofCluster->PosYMaxScal(0.75); //in % of length
+ Int_t iBRef = iCalSet % 1000;
+ Int_t iSet = (iCalSet - iBRef) / 1000;
+ Int_t iRSel = 0;
+ Int_t iRSelTyp = 0;
+ Int_t iRSelSm = 0;
+ Int_t iRSelRpc = 0;
+ iRSel = iBRef; // use diamond
+ Int_t iRSelin = iRSel;
+ iRSelRpc = iRSel % 10;
+ iRSelTyp = (iRSel - iRSelRpc) / 10;
+ iRSelSm = iRSelTyp % 10;
+ iRSelTyp = (iRSelTyp - iRSelSm) / 10;
+ tofCluster->SetBeamRefId(iRSelTyp); // define Beam reference counter
+ tofCluster->SetBeamRefSm(iRSelSm);
+ tofCluster->SetBeamRefDet(iRSelRpc);
+ tofCluster->SetBeamAddRefMul(-1);
+ tofCluster->SetBeamRefMulMax(3);
+ Int_t iSel2in = iSel2;
+ Int_t iSel2Rpc = iSel2 % 10;
+ iSel2 = (iSel2 - iSel2Rpc) / 10;
+ Int_t iSel2Sm = iSel2 % 10;
+ iSel2 = (iSel2 - iSel2Sm) / 10;
+ if (iSel2 > -1) {
+ tofCluster->SetSel2Id(iSel2);
+ tofCluster->SetSel2Sm(iSel2Sm);
+ tofCluster->SetSel2Rpc(iSel2Rpc);
+ }
+ Int_t iRef = iSet % 1000;
+ Int_t iDut = (iSet - iRef) / 1000;
+ Int_t iDutRpc = iDut % 10;
+ iDut = (iDut - iDutRpc) / 10;
+ Int_t iDutSm = iDut % 10;
+ iDut = (iDut - iDutSm) / 10;
+ tofCluster->SetDutId(iDut);
+ tofCluster->SetDutSm(iDutSm);
+ tofCluster->SetDutRpc(iDutRpc);
+
+ Int_t iRefRpc = iRef % 10;
+ iRef = (iRef - iRefRpc) / 10;
+ Int_t iRefSm = iRef % 10;
+ iRef = (iRef - iRefSm) / 10;
+
+ tofCluster->SetSelId(iRef);
+ tofCluster->SetSelSm(iRefSm);
+ tofCluster->SetSelRpc(iRefRpc);
+
+ run->AddTask(tofCluster);
+ std::cout << "-I- " << myName << ": Added task " << tofCluster->GetName() << std::endl;
+ } break;
+
+ default: {
+ ;
+ }
+ }
+ // -------------------------------------------------------------------------
+
+
+ // =========================================================================
+ // === Tof Tracking ===
+ // =========================================================================
+
+ cout << "<I> Initialize Tof tracker by ini_trks" << endl;
+ TString cTrkFile = Form("%s/%s_tofFindTracks.hst.root", TofFileFolder.Data(), cCalId.Data());
+
+ // ----- Local selection variables --------------------------------------
+
+ Int_t iRef = iSel % 1000;
+ Int_t iDut = (iSel - iRef) / 1000;
+ Int_t iDutRpc = iDut % 10;
+ iDut = (iDut - iDutRpc) / 10;
+ Int_t iDutSm = iDut % 10;
+ iDut = (iDut - iDutSm) / 10;
+ Int_t iBucRpc = 0;
+
+ // =========================================================================
+ // === Tracking ===
+ // =========================================================================
+
+ CbmTofTrackFinder* tofTrackFinder = new CbmTofTrackFinderNN();
+ tofTrackFinder->SetMaxTofTimeDifference(0.2); // in ns/cm
+ Int_t TrackerPar = 0;
+ switch (TrackerPar) {
+ case 0: // for full mTof setup
+ tofTrackFinder->SetTxLIM(0.3); // max slope dx/dz
+ tofTrackFinder->SetTyLIM(0.3); // max dev from mean slope dy/dz
+ tofTrackFinder->SetTxMean(0.); // mean slope dy/dz
+ tofTrackFinder->SetTyMean(0.); // mean slope dy/dz
+ break;
+ case 1: // for double stack test counters
+ tofTrackFinder->SetTxMean(0.21); // mean slope dy/dz
+ tofTrackFinder->SetTyMean(0.18); // mean slope dy/dz
+ tofTrackFinder->SetTxLIM(0.15); // max slope dx/dz
+ tofTrackFinder->SetTyLIM(0.18); // max dev from mean slope dy/dz
+ break;
+ }
+
+ CbmTofTrackFitter* tofTrackFitter = new CbmTofTrackFitterKF(0, 211);
+ TFitter* MyFit = new TFitter(1); // initialize Minuit
+ tofTrackFinder->SetFitter(tofTrackFitter);
+ CbmTofFindTracks* tofFindTracks = new CbmTofFindTracks("TOF Track Finder");
+ tofFindTracks->UseFinder(tofTrackFinder);
+ tofFindTracks->UseFitter(tofTrackFitter);
+ tofFindTracks->SetCalOpt(iCalOpt); // 1 - update offsets, 2 - update walk, 0 - bypass
+ tofFindTracks->SetCorMode(iGenCor); // valid options: 0,1,2,3,4,5,6, 10 - 19
+ tofFindTracks->SetTtTarg(0.042); // target value Mar2021, after T0 fix (double stack run 1058)
+ tofFindTracks->SetCalParFileName(cTrkFile); // Tracker parameter value file name
+ tofFindTracks->SetBeamCounter(5, 0, 0); // default beam counter
+ tofFindTracks->SetR0Lim(20.);
+ tofFindTracks->SetStationMaxHMul(30); // Max Hit Multiplicity in any used station
+
+ tofFindTracks->SetT0MAX(dScalFac); // in ns
+ tofFindTracks->SetSIGT(0.08); // default in ns
+ tofFindTracks->SetSIGX(0.3); // default in cm
+ tofFindTracks->SetSIGY(0.45); // default in cm
+ tofFindTracks->SetSIGZ(0.05); // default in cm
+ tofFindTracks->SetUseSigCalib(bUseSigCalib); // ignore resolutions in CalPar file
+ tofTrackFinder->SetSIGLIM(dChi2Lim2 * 2.); // matching window in multiples of chi2
+ tofTrackFinder->SetChiMaxAccept(dChi2Lim2); // max tracklet chi2
+
+ Int_t iMinNofHits = -1;
+ Int_t iNStations = 0;
+ Int_t iNReqStations = 3;
+
+ switch (iTrackingSetup) {
+ case 0: // bypass mode
+ iMinNofHits = -1;
+ iNStations = 1;
+ tofFindTracks->SetStation(0, 5, 0, 0); // Diamond
+ break;
+
+ case 1: // for calibration mode of full setup
+ {
+ Double_t dTsig = dScalFac * 0.03;
+ tofFindTracks->SetSIGT(dTsig); // allow for variable deviations in ns
+ }
+ iMinNofHits = 3;
+ iNStations = 32;
+ iNReqStations = 4;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(1, 0, 2, 2);
+ tofFindTracks->SetStation(2, 0, 1, 2);
+ tofFindTracks->SetStation(3, 0, 0, 2);
+ tofFindTracks->SetStation(4, 0, 2, 1);
+ tofFindTracks->SetStation(5, 0, 1, 1);
+ tofFindTracks->SetStation(6, 0, 0, 1);
+ tofFindTracks->SetStation(7, 0, 2, 3);
+ tofFindTracks->SetStation(8, 0, 1, 3);
+ tofFindTracks->SetStation(9, 0, 0, 3);
+ tofFindTracks->SetStation(10, 0, 2, 0);
+ tofFindTracks->SetStation(11, 0, 1, 0);
+ tofFindTracks->SetStation(12, 0, 0, 0);
+ tofFindTracks->SetStation(13, 0, 2, 4);
+ tofFindTracks->SetStation(14, 0, 1, 4);
+ tofFindTracks->SetStation(15, 0, 0, 4);
+ tofFindTracks->SetStation(16, 0, 4, 0);
+ tofFindTracks->SetStation(17, 0, 3, 0);
+ tofFindTracks->SetStation(18, 0, 4, 1);
+ tofFindTracks->SetStation(19, 0, 3, 1);
+ tofFindTracks->SetStation(20, 0, 4, 2);
+ tofFindTracks->SetStation(21, 0, 3, 2);
+ tofFindTracks->SetStation(22, 0, 4, 3);
+ tofFindTracks->SetStation(23, 0, 3, 3);
+ tofFindTracks->SetStation(24, 0, 4, 4);
+ tofFindTracks->SetStation(25, 0, 3, 4);
+ tofFindTracks->SetStation(26, 9, 0, 0);
+ tofFindTracks->SetStation(27, 9, 1, 0);
+ tofFindTracks->SetStation(28, 9, 0, 1);
+ tofFindTracks->SetStation(29, 9, 1, 1);
+ tofFindTracks->SetStation(30, 6, 0, 0);
+ tofFindTracks->SetStation(31, 6, 0, 1);
+ break;
+
+ case 11: // for calibration mode of 2-stack & test counters
+ iMinNofHits = 4;
+ iNStations = 9;
+ iNReqStations = 5;
+ tofFindTracks->SetStation(0, 0, 4, 1);
+ tofFindTracks->SetStation(1, 9, 0, 0);
+ tofFindTracks->SetStation(2, 9, 1, 0);
+ tofFindTracks->SetStation(3, 9, 0, 1);
+ tofFindTracks->SetStation(4, 9, 1, 1);
+ tofFindTracks->SetStation(5, 0, 3, 1);
+ tofFindTracks->SetStation(6, 0, 4, 0);
+ tofFindTracks->SetStation(7, 0, 3, 2);
+ tofFindTracks->SetStation(8, 5, 0, 0);
+ break;
+
+ case 2:
+ iMinNofHits = 5;
+ iNStations = 28;
+ iNReqStations = 5;
+ tofFindTracks->SetStation(0, 0, 2, 2);
+ tofFindTracks->SetStation(1, 0, 0, 2);
+ tofFindTracks->SetStation(2, 0, 1, 2);
+ tofFindTracks->SetStation(3, 0, 2, 1);
+ tofFindTracks->SetStation(4, 0, 0, 1);
+ tofFindTracks->SetStation(5, 0, 1, 1);
+ tofFindTracks->SetStation(6, 0, 2, 3);
+ tofFindTracks->SetStation(7, 0, 0, 3);
+ tofFindTracks->SetStation(8, 0, 1, 3);
+ tofFindTracks->SetStation(9, 0, 2, 0);
+ tofFindTracks->SetStation(10, 0, 0, 0);
+ tofFindTracks->SetStation(11, 0, 1, 0);
+ tofFindTracks->SetStation(12, 0, 2, 4);
+ tofFindTracks->SetStation(13, 0, 0, 4);
+ tofFindTracks->SetStation(14, 0, 1, 4);
+ tofFindTracks->SetStation(15, 0, 4, 0);
+ tofFindTracks->SetStation(16, 0, 3, 0);
+ tofFindTracks->SetStation(17, 0, 4, 1);
+ tofFindTracks->SetStation(18, 0, 3, 1);
+ tofFindTracks->SetStation(19, 0, 4, 2);
+ tofFindTracks->SetStation(20, 0, 3, 2);
+ tofFindTracks->SetStation(21, 0, 4, 3);
+ tofFindTracks->SetStation(22, 0, 3, 3);
+ tofFindTracks->SetStation(23, 0, 4, 4);
+ tofFindTracks->SetStation(24, 0, 3, 4);
+ tofFindTracks->SetStation(25, 9, 0, 0);
+ tofFindTracks->SetStation(26, 9, 0, 1);
+ tofFindTracks->SetStation(27, 5, 0, 0);
+ break;
+
+ case 3:
+ iMinNofHits = 3;
+ iNStations = 16;
+ iNReqStations = 4;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(1, 0, 2, 2);
+ tofFindTracks->SetStation(2, 0, 1, 2);
+ tofFindTracks->SetStation(3, 0, 0, 2);
+
+ tofFindTracks->SetStation(4, 0, 2, 1);
+ tofFindTracks->SetStation(5, 0, 1, 1);
+ tofFindTracks->SetStation(6, 0, 0, 1);
+
+ tofFindTracks->SetStation(7, 0, 2, 3);
+ tofFindTracks->SetStation(8, 0, 1, 3);
+ tofFindTracks->SetStation(9, 0, 0, 3);
+
+ tofFindTracks->SetStation(10, 0, 2, 0);
+ tofFindTracks->SetStation(11, 0, 1, 0);
+ tofFindTracks->SetStation(12, 0, 0, 0);
+
+ tofFindTracks->SetStation(13, 0, 2, 4);
+ tofFindTracks->SetStation(14, 0, 1, 4);
+ tofFindTracks->SetStation(15, 0, 0, 4);
+
+ /*
+ tofFindTracks->SetStation(16, 0, 3, 2);
+ tofFindTracks->SetStation(17, 0, 4, 2);
+ tofFindTracks->SetStation(18, 0, 3, 1);
+ tofFindTracks->SetStation(19, 0, 4, 1);
+ tofFindTracks->SetStation(20, 0, 3, 3);
+ tofFindTracks->SetStation(21, 0, 4, 3);
+ tofFindTracks->SetStation(22, 0, 3, 0);
+ tofFindTracks->SetStation(23, 0, 4, 0);
+ tofFindTracks->SetStation(24, 0, 3, 4);
+ tofFindTracks->SetStation(25, 0, 4, 4);
+ */
+ break;
+
+ case 4: // for USTC evaluation (dut=910,911)
+ iMinNofHits = 4;
+ iNStations = 6;
+ iNReqStations = 6;
+ tofFindTracks->SetStation(0, 0, 4, 1);
+ tofFindTracks->SetStation(1, 0, 3, 1);
+ tofFindTracks->SetStation(2, 9, 0, 1);
+ tofFindTracks->SetStation(3, 9, 0, 0);
+ tofFindTracks->SetStation(4, 5, 0, 0);
+ tofFindTracks->SetStation(5, iDut, iDutSm, iDutRpc);
+ break;
+
+ case 14:
+ iMinNofHits = 3;
+ iNStations = 15;
+ iNReqStations = 4;
+ tofFindTracks->SetStation(0, 0, 2, 2);
+ tofFindTracks->SetStation(1, 0, 1, 2);
+ tofFindTracks->SetStation(2, 0, 0, 2);
+ tofFindTracks->SetStation(0, 0, 2, 1);
+ tofFindTracks->SetStation(1, 0, 1, 1);
+ tofFindTracks->SetStation(2, 0, 0, 1);
+ tofFindTracks->SetStation(0, 0, 2, 0);
+ tofFindTracks->SetStation(1, 0, 1, 0);
+ tofFindTracks->SetStation(2, 0, 0, 0);
+ tofFindTracks->SetStation(0, 0, 2, 3);
+ tofFindTracks->SetStation(1, 0, 1, 3);
+ tofFindTracks->SetStation(2, 0, 0, 3);
+ tofFindTracks->SetStation(0, 0, 2, 4);
+ tofFindTracks->SetStation(1, 0, 1, 4);
+ tofFindTracks->SetStation(2, 0, 0, 4);
+ break;
+
+ case 5: // for calibration of 2-stack and add-on counters (STAR2, BUC)
+ iMinNofHits = 3;
+ iNStations = 7;
+ iNReqStations = 4;
+ tofFindTracks->SetStation(6, 0, 4, 1);
+ tofFindTracks->SetStation(1, 6, 0, 1);
+ tofFindTracks->SetStation(2, 9, 0, 0);
+ tofFindTracks->SetStation(3, 9, 0, 1);
+ tofFindTracks->SetStation(4, 6, 0, 0);
+ tofFindTracks->SetStation(5, 0, 3, 1);
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ break;
+
+ case 6: // for double stack USTC counter evaluation
+ iMinNofHits = 5;
+ iNStations = 6;
+ iNReqStations = 6;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(1, 6, 0, 1);
+ tofFindTracks->SetStation(2, 0, 4, 1);
+ tofFindTracks->SetStation(3, 6, 0, 0);
+ tofFindTracks->SetStation(4, 0, 3, 1);
+ tofFindTracks->SetStation(5, iDut, iDutSm, iDutRpc);
+ break;
+
+ case 7: // for double stack USTC counter evaluation
+ iMinNofHits = 3;
+ iNStations = 4;
+ iNReqStations = 4;
+ tofFindTracks->SetStation(0, 0, 4, 1);
+ tofFindTracks->SetStation(1, 6, 0, 1);
+ tofFindTracks->SetStation(2, 6, 0, 0);
+ tofFindTracks->SetStation(3, iDut, iDutSm, iDutRpc);
+ break;
+
+ case 8: // evaluation of add-on counters (BUC)
+ iMinNofHits = 5;
+ iNStations = 6;
+ iNReqStations = 6;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(1, 9, 0, 1);
+ tofFindTracks->SetStation(2, 0, 4, 1);
+ tofFindTracks->SetStation(3, 9, 0, 0);
+ tofFindTracks->SetStation(4, 0, 3, 1);
+ tofFindTracks->SetStation(5, iDut, iDutSm, iDutRpc);
+ break;
+
+ case 9: // calibration of Star2
+ iMinNofHits = 4;
+ iNStations = 5;
+ iNReqStations = 5;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(2, 9, 0, 1);
+ tofFindTracks->SetStation(1, 0, 4, 1);
+ tofFindTracks->SetStation(3, 9, 0, 0);
+ tofFindTracks->SetStation(4, 0, 3, 1);
+ break;
+
+ case 10:
+ iMinNofHits = 3;
+ iNStations = 4;
+ iNReqStations = 4;
+ tofFindTracks->SetStation(0, 5, 0, 0);
+ tofFindTracks->SetStation(3, 0, 1, 2);
+ tofFindTracks->SetStation(2, 0, 0, 2);
+ tofFindTracks->SetStation(1, 0, 2, 2);
+ break;
+
+ default:
+ cout << "Tracking setup " << iTrackingSetup << " not implemented " << endl;
+ return 1;
+ ;
+ }
+ tofFindTracks->SetMinNofHits(iMinNofHits);
+ tofFindTracks->SetNStations(iNStations);
+ tofFindTracks->SetNReqStations(iNReqStations);
+ tofFindTracks->PrintSetup();
+ std::cout << "MinNofHitsPerTrack: " << iMinNofHits << std::endl;
+ run->AddTask(tofFindTracks);
+ }
+
+
+ // =========================================================================
+ // === L1 ===
+ // =========================================================================
+
+ // CbmKF* kalman = new CbmKF();
+ // run->AddTask(kalman);
+ // CbmL1* l1 = new CbmL1();
+ // l1->SetLegacyEventMode(1);
+ // l1->SetMcbmMode();
+ // l1->SetUseHitErrors(1);
+ // if (strcmp(geoSetupTag.data(), "mcbm_beam_2021_07_surveyed") == 0) l1->SetMissingHits(1);
+ //
+ // // --- Material budget file names
+ // TString mvdGeoTag;
+ // if (geoSetup->GetGeoTag(ECbmModuleId::kMvd, mvdGeoTag)) {
+ // TString parFile = gSystem->Getenv("VMCWORKDIR");
+ // parFile = parFile + "/parameters/mvd/mvd_matbudget_" + mvdGeoTag + ".root";
+ // std::cout << "Using material budget file " << parFile << std::endl;
+ // l1->SetMvdMaterialBudgetFileName(parFile.Data());
+ // }
+ // TString stsGeoTag;
+ // if (geoSetup->GetGeoTag(ECbmModuleId::kSts, stsGeoTag)) {
+ // TString parFile = gSystem->Getenv("VMCWORKDIR");
+ // parFile = parFile + "/parameters/sts/sts_matbudget_v19a.root";
+ // std::cout << "Using material budget file " << parFile << std::endl;
+ // l1->SetStsMaterialBudgetFileName(parFile.Data());
+ // }
+ //
+ // 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 parFile = gSystem->Getenv("VMCWORKDIR");
+ // parFile = parFile + "/parameters/much/much_" + geoTag(0, 4) + "_digi_sector.root";
+ // std::cout << "L1: Using parameter file " << parFile << std::endl;
+ // l1->SetMuchPar(parFile);
+ //
+ // TString parFile2 = gSystem->Getenv("VMCWORKDIR");
+ // parFile2 = parFile2 + "/parameters/much/much_matbudget_" + geoTag + ".root ";
+ // std::cout << "Using material budget file " << parFile2 << std::endl;
+ // l1->SetMuchMaterialBudgetFileName(parFile2.Data());
+ // }
+ //
+ // TString trdGeoTag;
+ // if (geoSetup->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());
+ // }
+ //
+ // TString tofGeoTag;
+ // if (geoSetup->GetGeoTag(ECbmModuleId::kTof, tofGeoTag)) {
+ // TString parFile = gSystem->Getenv("VMCWORKDIR");
+ // parFile = parFile + "/parameters/tof/tof_matbudget_" + tofGeoTag + ".root ";
+ // std::cout << "Using material budget file " << parFile << std::endl;
+ // l1->SetTofMaterialBudgetFileName(parFile.Data());
+ // }
+ //
+
+ // Workaround to get it running:
+ // 1) Change fUseGlobal in line 129 of CbmStsParSetModule.h to
+ // Bool_t fUseGlobal = kTRUE;
+ // 2) Change fUseGlobal in line 114 of CbmStsParSetSensor.h to
+ // Bool_t fUseGlobal = kTRUE;
+ //run->AddTask(l1);
+
+ // CbmL1GlobalTrackFinder* globalTrackFinder = new CbmL1GlobalTrackFinder();
+ // FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder);
+ //run->AddTask(globalFindTracks);
+
+
+ // ----- 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);
+ rtdb->saveOutput();
+ rtdb->print();
+
+ // ----- Run initialisation -------------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+
+ // ------------------------------------------------------------------------
+
+ // =========================================================================
+ // === QA ===
+ // =========================================================================
+
+ // e.g for RICH:
+ // CbmRichMCbmQaReal* qaTask = new CbmRichMCbmQaReal();
+ // if (taskId < 0) { qaTask->SetOutputDir(Form("result_run%d", uRunId)); }
+ // else {
+ // qaTask->SetOutputDir(Form("result_run%d_%05d", uRunId, taskId));
+ // }
+ // qaTask->XOffsetHistos(+25.3);
+ // qaTask->SetMaxNofDrawnEvents(100);
+ // qaTask->SetTotRich(23.7, 30.0);
+ // qaTask->SetTriggerRichHits(eb_TriggerMinNumberRich);
+ // qaTask->SetTriggerTofHits(1);
+ // run->AddTask(qaTask);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Start run ----------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, 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 " << 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;
+}
diff --git a/macro/beamtime/mcbm2021/mcbm_inspect_alignment.C b/macro/beamtime/mcbm2021/mcbm_inspect_alignment.C
new file mode 100644
index 0000000000..d05efabe7f
--- /dev/null
+++ b/macro/beamtime/mcbm2021/mcbm_inspect_alignment.C
@@ -0,0 +1,1214 @@
+/* Copyright (C) 2022 UGiessen, Giessen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Adrian Weber [committer]*/
+
+void initHists(CbmHistManager* fHM, Int_t nofEvents);
+void drawHists(CbmHistManager* fHM);
+TVector3 extrpolate(CbmTofHit* tofHit, TVector3* vertex, Double_t Z);
+
+void ProjectXY(std::string name, CbmHistManager* fHM, std::vector<TH1*>* HistoX, std::vector<std::string>* stringX,
+ std::vector<TH1*>* HistoY, std::vector<std::string>* stringY);
+
+void ProjectAxis(std::string name, CbmHistManager* fHM, char axis, std::vector<TH1*>* Histo,
+ std::vector<std::string>* string);
+
+std::string FitTimeOffset(TH1* hist);
+
+struct ClosestCbmHit {
+ Bool_t Existing = false;
+ Double_t X = 10000.0;
+ Double_t Y = 10000.0;
+ Double_t Z = 10000.0;
+ Double_t Dist = 10000.0;
+ size_t Indx = -1;
+ Double_t X_pos = 10000.0;
+ Double_t Y_pos = 10000.0;
+};
+
+char outdir[256];
+
+void mcbm_inspect_alignment(Int_t runId = 1588, Int_t nofEvents = 200)
+{
+
+ CbmHistManager* fHM = new CbmHistManager();
+
+ TTree* InputTreeRECO;
+ TFile* fRECO;
+
+ TVector3 vertexPoint(-0.0, 0.0, 0.0);
+ Double_t CutCloseDistStsTrack = 60.0; // cm
+
+ char name[256], dir[256];
+
+ //sprintf(dir, "/data/cbmroot/cbmsource/macro/beamtime/mcbm2021/rec/digi5");
+ //sprintf(dir, "/data/cbmroot/cbmsource/macro/beamtime/mcbm2021/rec/aligndigi5");
+ //sprintf(dir, "/data/cbmroot/cbmsource/macro/beamtime/mcbm2021/rec/aligndigi5_2HitTrack");
+ //sprintf(dir, "/data/cbmroot/cbmsource/macro/beamtime/mcbm2021/rec/aligndigi2_2HitTrack");
+ //sprintf(dir, "/data/cbmroot/cbmsource/macro/beamtime/mcbm2021/rec/tofPar0");
+ sprintf(dir, "/data/cbmroot/cbmsource/macro/beamtime/mcbm2021/rec/tofPar2_noAlign_L1");
+ //sprintf(dir, "/data/cbmroot/cbmsource/macro/beamtime/mcbm2021/rec/tofPar0_noAlign_5hits");
+
+ sprintf(name, "%s/reco_event_mcbm_%u.root", dir, runId);
+
+ sprintf(outdir, "./results_mTofPar2_AlignmentMacro");
+
+ std::cout << name << std::endl;
+
+ initHists(fHM, nofEvents);
+
+ auto tofHitYMin = 5000.0;
+ auto tofHitYMax = -5000.0;
+
+ fRECO = new TFile(name);
+ if (!fRECO || fRECO->IsZombie() || fRECO->GetNkeys() < 1 || fRECO->TestBit(TFile::kRecovered)) { fRECO->Close(); }
+
+ InputTreeRECO = (TTree*) fRECO->Get("cbmsim");
+ if (InputTreeRECO == nullptr) std::cout << "No cbmsim found" << std::endl;
+
+ TClonesArray* fTofHit = nullptr;
+ TClonesArray* fTofTrack = nullptr;
+ TClonesArray* fStsHit = nullptr;
+ TClonesArray* fCbmEvent = nullptr;
+ InputTreeRECO->SetBranchAddress("TofHit", &fTofHit);
+ InputTreeRECO->SetBranchAddress("TofTracks", &fTofTrack);
+ InputTreeRECO->SetBranchAddress("StsHit", &fStsHit);
+ InputTreeRECO->SetBranchAddress("CbmEvent", &fCbmEvent);
+
+ auto nTS = InputTreeRECO->GetEntries();
+ if (nofEvents > 0 && nTS >= nofEvents) { nTS = nofEvents; }
+ std::cout << "Entries: " << InputTreeRECO->GetEntries() << std::endl;
+
+ TFile* oldFileLin = gFile;
+ TDirectory* oldirLin = gDirectory;
+
+ std::string s = Form("%s/GraphFitsLinear.root", outdir);
+ TFile* outFileLin = new TFile(s.c_str(), "RECREATE");
+
+
+ for (int iTS = 0; iTS < nTS; iTS++) {
+ fHM->H1("fhNofTimeslices")->Fill(1);
+ InputTreeRECO->GetEntry(iTS);
+ std::cout << "Timeslice No. " << iTS << std::endl;
+ Int_t nTParts = fTofHit->GetEntriesFast();
+ Int_t nEv = fCbmEvent->GetEntriesFast();
+
+ std::cout << "Events in TS: " << nEv << std::endl;
+ fHM->H1("fhEventsPerTS")->Fill(iTS, nEv);
+
+ for (auto iEv = 0; iEv < nEv; iEv++) {
+ CbmEvent* ev = static_cast<CbmEvent*>(fCbmEvent->At(iEv));
+ auto eventStartTime = ev->GetStartTime();
+
+ // Hit-Multiplicty Cuts:
+ if (ev->GetNofData(ECbmDataType::kTofHit) < 2) continue;
+
+ // Fill Hit Multiplicity:
+ fHM->H1("fhMultiplicityStsHits")->Fill(ev->GetNofData(ECbmDataType::kStsHit));
+ fHM->H1("fhMultiplicityTofHits")->Fill(ev->GetNofData(ECbmDataType::kTofHit));
+
+
+ //Correct the TofHits
+ std::vector<size_t> tofIndxStation[2];
+ for (int j = 0; j < ev->GetNofData(ECbmDataType::kTofHit); j++) {
+ auto iTofHit = ev->GetIndex(ECbmDataType::kTofHit, j);
+ CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHit->At(iTofHit));
+ auto TofStationNo = (tofHit->GetAddress() >> 4) & 0xf;
+
+ if (TofStationNo == 0) { tofIndxStation[0].push_back(iTofHit); }
+ else {
+ tofIndxStation[1].push_back(iTofHit);
+ }
+ } // TofHits
+
+ for (int j = 0; j < tofIndxStation[1].size(); j++) {
+ CbmTofHit* tofHit_stat1 = static_cast<CbmTofHit*>(fTofHit->At(tofIndxStation[1].at(j)));
+
+ ClosestCbmHit closestTofHits; // Station 0
+ for (int k = 0; k < tofIndxStation[0].size(); k++) {
+ CbmTofHit* tofHit_stat0 = static_cast<CbmTofHit*>(fTofHit->At(tofIndxStation[0].at(k)));
+ TVector3 vExtr = extrpolate(tofHit_stat1, &vertexPoint, tofHit_stat0->GetZ());
+
+ auto Xextr = vExtr.X(); //tofHit->GetX()* stsHit->GetZ() / tofHit->GetZ();
+ auto Yextr = vExtr.Y(); //tofHit->GetY()* stsHit->GetZ() / tofHit->GetZ();
+
+ auto Xdiff = Xextr - tofHit_stat0->GetX();
+ auto Ydiff = Yextr - tofHit_stat0->GetY();
+
+ Double_t distTof = TMath::Sqrt(Xdiff * Xdiff + Ydiff * Ydiff);
+
+ if (distTof < closestTofHits.Dist && 0.5 > distTof) {
+ closestTofHits.Existing = true;
+ closestTofHits.X = Xdiff;
+ closestTofHits.Y = Ydiff;
+ closestTofHits.Z = tofHit_stat0->GetZ();
+ closestTofHits.Dist = distTof;
+ closestTofHits.Indx = tofIndxStation[0].at(k);
+ closestTofHits.X_pos = tofHit_stat0->GetX();
+ closestTofHits.Y_pos = tofHit_stat0->GetY();
+ }
+ } // TofHits
+
+ ClosestCbmHit closestStsHitArray[3];
+ ClosestCbmHit closestStsHitStations[2];
+
+ auto nofStsHitsEv = ev->GetNofData(ECbmDataType::kStsHit);
+ for (int jSts = 0; jSts < nofStsHitsEv; jSts++) {
+ auto iStsHit = ev->GetIndex(ECbmDataType::kStsHit, jSts);
+ CbmStsHit* stsHit = static_cast<CbmStsHit*>(fStsHit->At(iStsHit));
+
+ TVector3 vExtr = extrpolate(tofHit_stat1, &vertexPoint, stsHit->GetZ());
+
+ auto Xextr = vExtr.X(); //tofHit->GetX()* stsHit->GetZ() / tofHit->GetZ();
+ auto Yextr = vExtr.Y(); //tofHit->GetY()* stsHit->GetZ() / tofHit->GetZ();
+
+ auto stsUnit = 0;
+ auto stsStation = (stsHit->GetAddress() >> 4) & 0xF;
+ auto stsUnitAdd = (stsHit->GetAddress() >> 8) & 0x4;
+
+ // if (stsHit->GetZ() < 28.0+8.5){
+ if (stsStation == 0 && stsUnitAdd == 0) { stsUnit = 0; }
+ else if (stsStation == 0 && stsUnitAdd != 0) {
+ stsUnit = 1;
+ }
+ else {
+ stsStation = 1;
+ stsUnit = 2;
+ }
+
+ auto Xdiff = Xextr - stsHit->GetX();
+ auto Ydiff = Yextr - stsHit->GetY();
+
+ Double_t distSts = TMath::Sqrt(Xdiff * Xdiff + Ydiff * Ydiff);
+
+ if (distSts < closestStsHitArray[stsUnit].Dist) {
+ closestStsHitArray[stsUnit].Existing = true;
+ closestStsHitArray[stsUnit].X = Xdiff;
+ closestStsHitArray[stsUnit].Y = Ydiff;
+ closestStsHitArray[stsUnit].Z = stsHit->GetZ();
+ closestStsHitArray[stsUnit].Dist = distSts;
+ closestStsHitArray[stsUnit].Indx = iStsHit;
+ closestStsHitArray[stsUnit].X_pos = stsHit->GetX();
+ closestStsHitArray[stsUnit].Y_pos = stsHit->GetY();
+ }
+
+ if (distSts < closestStsHitStations[stsStation].Dist) {
+ closestStsHitStations[stsStation].Existing = true;
+ closestStsHitStations[stsStation].X = Xdiff;
+ closestStsHitStations[stsStation].Y = Ydiff;
+ closestStsHitStations[stsStation].Z = stsHit->GetZ();
+ closestStsHitStations[stsStation].Dist = distSts;
+ closestStsHitStations[stsStation].Indx = iStsHit;
+ closestStsHitStations[stsStation].X_pos = stsHit->GetX();
+ closestStsHitStations[stsStation].Y_pos = stsHit->GetY();
+ }
+ } //STS Loop
+
+ std::vector<Double_t> ZArray, XArray, YArray, IndexArray;
+ ZArray.push_back(tofHit_stat1->GetZ());
+ XArray.push_back(tofHit_stat1->GetX());
+ YArray.push_back(tofHit_stat1->GetY());
+ IndexArray.push_back(tofIndxStation[1].at(j));
+ if (closestTofHits.Existing) {
+ ZArray.push_back(closestTofHits.Z);
+ XArray.push_back(closestTofHits.X_pos);
+ YArray.push_back(closestTofHits.Y_pos);
+ IndexArray.push_back(closestTofHits.Indx);
+ }
+ if (closestStsHitStations[0].Existing) {
+ ZArray.push_back(closestStsHitStations[0].Z);
+ XArray.push_back(closestStsHitStations[0].X_pos);
+ YArray.push_back(closestStsHitStations[0].Y_pos);
+ IndexArray.push_back(closestStsHitStations[0].Indx);
+ }
+ if (closestStsHitStations[1].Existing) {
+ ZArray.push_back(closestStsHitStations[1].Z);
+ XArray.push_back(closestStsHitStations[1].X_pos);
+ YArray.push_back(closestStsHitStations[1].Y_pos);
+ IndexArray.push_back(closestStsHitStations[1].Indx);
+ }
+
+ if (ZArray.size() > 3) {
+ //std::cout << "Event # " << iEv << std::endl;
+ //std::cout << "Xpos: "<< closestTofHits.X_pos<< std::endl;
+ TGraph* grX = new TGraph(ZArray.size(), &ZArray[0], &XArray[0]);
+ TF1* linFitX = new TF1("f1", "[0]+[1]*x", 0, 0.1);
+ grX->Fit(linFitX, "Q");
+
+ TGraph* grY = new TGraph(ZArray.size(), &ZArray[0], &YArray[0]);
+ TF1* linFitY = new TF1("f2", "[0]+[1]*x", 0, 0.1);
+ grY->Fit(linFitY, "Q");
+
+ fHM->H1("fhLinearFitChi2X")->Fill(linFitX->GetChisquare());
+ fHM->H1("fhLinearFitChi2Y")->Fill(linFitY->GetChisquare());
+
+ if (linFitX->GetChisquare() < 0.5 && linFitY->GetChisquare() < 0.5) {
+ fHM->H2("fhLinearFitXY")->Fill(linFitX->GetParameter(0), linFitY->GetParameter(0));
+ fHM->H1("fhLinearFitZX")->Fill(linFitX->GetParameter(0));
+ fHM->H1("fhLinearFitZY")->Fill(linFitY->GetParameter(0));
+
+ fHM->H2("fhLinearFitXY_.25")
+ ->Fill(linFitX->GetParameter(0) + (.25) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (.25) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_.50")
+ ->Fill(linFitX->GetParameter(0) + (.50) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (.50) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_.75")
+ ->Fill(linFitX->GetParameter(0) + (.75) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (.75) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_1")
+ ->Fill(linFitX->GetParameter(0) + (1) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (1) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_2")
+ ->Fill(linFitX->GetParameter(0) + (2) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (2) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_4")
+ ->Fill(linFitX->GetParameter(0) + (4) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (4) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_6")
+ ->Fill(linFitX->GetParameter(0) + (6) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (6) * linFitY->GetParameter(1));
+
+ fHM->H2("fhLinearFitXY_-.5")
+ ->Fill(linFitX->GetParameter(0) + (-.5) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (-.5) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_-1")
+ ->Fill(linFitX->GetParameter(0) + (-1) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (-1) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_-2")
+ ->Fill(linFitX->GetParameter(0) + (-2) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (-2) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_-4")
+ ->Fill(linFitX->GetParameter(0) + (-4) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (-4) * linFitY->GetParameter(1));
+ fHM->H2("fhLinearFitXY_-6")
+ ->Fill(linFitX->GetParameter(0) + (-6) * linFitX->GetParameter(1),
+ linFitY->GetParameter(0) + (-6) * linFitY->GetParameter(1));
+
+ // Good Fit!: Extract informations of this kind of tracks
+ const CbmTofHit* tofHit_sel_0 = static_cast<CbmTofHit*>(fTofHit->At(IndexArray.at(0)));
+ const CbmTofHit* tofHit_sel_1 = static_cast<CbmTofHit*>(fTofHit->At(IndexArray.at(1)));
+
+ auto selectedTrackTime = tofHit_sel_0->GetTime() - tofHit_sel_1->GetTime();
+ auto selectedTrackXdiff = tofHit_sel_0->GetX() - tofHit_sel_1->GetX();
+ auto selectedTrackYdiff = tofHit_sel_0->GetY() - tofHit_sel_1->GetY();
+ auto selectedTrackZdiff = tofHit_sel_0->GetZ() - tofHit_sel_1->GetZ();
+ auto Rdiff = TMath::Sqrt(selectedTrackXdiff * selectedTrackXdiff + selectedTrackYdiff * selectedTrackYdiff
+ + selectedTrackZdiff * selectedTrackZdiff);
+ auto selectedTrackBeta = (1e7 / (selectedTrackTime / Rdiff)) / TMath::C();
+
+ fHM->H1("fhSimpleTrackTofTime")->Fill(selectedTrackTime);
+ fHM->H1("fhSimpleTrackTofDist")->Fill(Rdiff);
+ fHM->H1("fhSimpleTrackTofBeta")->Fill(selectedTrackBeta);
+ }
+
+ //linFitX->Print();
+ if (outFileLin->IsOpen()) {
+ if (iEv < 1000 && linFitX->GetChisquare() < .5) grX->Write();
+ }
+ }
+ } // TofHits
+
+
+ std::vector<size_t> stsIndxStation[2];
+ for (int jSts = 0; jSts < ev->GetNofData(ECbmDataType::kStsHit); jSts++) {
+
+ auto iStsHit = ev->GetIndex(ECbmDataType::kStsHit, jSts);
+ CbmStsHit* stsHit = static_cast<CbmStsHit*>(fStsHit->At(iStsHit));
+ auto x = stsHit->GetX();
+ auto y = stsHit->GetY();
+ auto z = stsHit->GetZ();
+
+ fHM->H1("fhStsHitsZ")->Fill(z);
+
+ auto stsStation = (stsHit->GetAddress() >> 4) & 0xF;
+ auto stsUnit = (stsHit->GetAddress() >> 8) & 0x4;
+
+ if (stsStation == 0 && stsUnit == 0) {
+ fHM->H2("fhStsUnit0Hits")->Fill(stsHit->GetX(), stsHit->GetY());
+ fHM->H2("fhStsStation0Hits")->Fill(stsHit->GetX(), stsHit->GetY());
+ stsIndxStation[0].push_back(iStsHit);
+ }
+ else if (stsStation == 0 && stsUnit != 0) {
+ fHM->H2("fhStsUnit1Hits")->Fill(stsHit->GetX(), stsHit->GetY());
+ fHM->H2("fhStsStation0Hits")->Fill(stsHit->GetX(), stsHit->GetY());
+ stsIndxStation[0].push_back(iStsHit);
+ }
+ else {
+ fHM->H2("fhStsUnit2Hits")->Fill(stsHit->GetX(), stsHit->GetY());
+ stsIndxStation[1].push_back(iStsHit);
+ }
+ } // StsHits
+
+
+ auto nofTofHitsEv = ev->GetNofData(ECbmDataType::kTofHit);
+ for (int j = 0; j < nofTofHitsEv; j++) {
+ auto iTofHit = ev->GetIndex(ECbmDataType::kTofHit, j);
+ CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHit->At(iTofHit));
+
+ fHM->H1("fhTofHitsTime")->Fill(tofHit->GetTime() - eventStartTime);
+
+ auto TofStationNo = (tofHit->GetAddress() >> 4) & 0xf;
+ auto modPos = (tofHit->GetAddress() >> 15) & 0x1; // 0 is closer to vertex
+ auto modPosY = (tofHit->GetAddress() >> 16) & 0xf; // ModPos 0: 0,1,2 ModPos 1: 0,1
+
+ fHM->H1("fhTofHitsZ")->Fill(tofHit->GetZ());
+ fHM->H2("fhTofHitsXY")->Fill(tofHit->GetX(), tofHit->GetY());
+
+
+ if (TofStationNo == 0 && modPos == 0 && modPosY == 1) {
+ fHM->H2("fhTofHitsXY_001")->Fill(tofHit->GetX(), tofHit->GetY());
+ fHM->H1("fhTofHitsTime_001")->Fill(tofHit->GetTime() - eventStartTime);
+ if (tofHitYMin > tofHit->GetY()) tofHitYMin = tofHit->GetY();
+ if (tofHitYMax < tofHit->GetY()) tofHitYMax = tofHit->GetY();
+ }
+
+ if (TofStationNo == 0) { fHM->H1("fhTofHitsTimeLayer0")->Fill(tofHit->GetTime() - eventStartTime); }
+
+ if (TofStationNo == 1) { fHM->H1("fhTofHitsTimeLayer1")->Fill(tofHit->GetTime() - eventStartTime); }
+
+ ClosestCbmHit closestStsHitArray[3];
+ ClosestCbmHit closestStsHitStations[2];
+
+ auto nofStsHitsEv = ev->GetNofData(ECbmDataType::kStsHit);
+ for (int jSts = 0; jSts < nofStsHitsEv; jSts++) {
+ auto iStsHit = ev->GetIndex(ECbmDataType::kStsHit, jSts);
+ CbmStsHit* stsHit = static_cast<CbmStsHit*>(fStsHit->At(iStsHit));
+
+ TVector3 vExtr = extrpolate(tofHit, &vertexPoint, stsHit->GetZ());
+
+ auto Xextr = vExtr.X(); //tofHit->GetX()* stsHit->GetZ() / tofHit->GetZ();
+ auto Yextr = vExtr.Y(); //tofHit->GetY()* stsHit->GetZ() / tofHit->GetZ();
+
+ auto stsUnit = 0;
+
+ auto stsStation = (stsHit->GetAddress() >> 4) & 0xF;
+ auto stsUnitAdd = (stsHit->GetAddress() >> 8) & 0x4;
+
+ if (stsStation == 0 && stsUnitAdd == 0) { stsUnit = 0; }
+ else if (stsStation == 0 && stsUnitAdd != 0) {
+ stsUnit = 1;
+ }
+ else {
+ stsStation = 1;
+ stsUnit = 2;
+ }
+
+ auto Xdiff = Xextr - stsHit->GetX();
+ auto Ydiff = Yextr - stsHit->GetY();
+
+ Double_t distSts = TMath::Sqrt(Xdiff * Xdiff + Ydiff * Ydiff);
+
+ if (distSts < closestStsHitArray[stsUnit].Dist) {
+ closestStsHitArray[stsUnit].Existing = true;
+ closestStsHitArray[stsUnit].X = Xdiff;
+ closestStsHitArray[stsUnit].Y = Ydiff;
+ closestStsHitArray[stsUnit].Z = stsHit->GetZ();
+ closestStsHitArray[stsUnit].Dist = distSts;
+ closestStsHitArray[stsUnit].Indx = iStsHit;
+ closestStsHitArray[stsUnit].X_pos = stsHit->GetX();
+ closestStsHitArray[stsUnit].Y_pos = stsHit->GetY();
+ }
+
+ if (distSts < closestStsHitStations[stsStation].Dist) {
+ closestStsHitStations[stsStation].Existing = true;
+ closestStsHitStations[stsStation].X = Xdiff;
+ closestStsHitStations[stsStation].Y = Ydiff;
+ closestStsHitStations[stsStation].Z = stsHit->GetZ();
+ closestStsHitStations[stsStation].Dist = distSts;
+ closestStsHitStations[stsStation].Indx = iStsHit;
+ closestStsHitStations[stsStation].X_pos = stsHit->GetX();
+ closestStsHitStations[stsStation].Y_pos = stsHit->GetY();
+ }
+
+
+ fHM->H2(Form("fhTofSts%uXYdiff", stsUnit))->Fill(Xdiff, Ydiff);
+
+ if (TMath::Abs(Xdiff) < 0.2 && TMath::Abs(Ydiff) < 0.2) {
+ vExtr = extrpolate(tofHit, &vertexPoint, 8.5);
+ fHM->H2(Form("fhTofSts%uXY_8.5cm", stsUnit))->Fill(vExtr.X(), vExtr.Y());
+ }
+ } //STS Loop
+
+ for (auto o = 0; o < 3; ++o) {
+ if (closestStsHitArray[o].Existing == true && (closestStsHitStations[0].Existing == true)
+ && (closestStsHitStations[1].Existing == true)) {
+ TVector3 vExtr_tof = extrpolate(tofHit, &vertexPoint, closestStsHitArray[o].Z);
+ fHM->H2(Form("fhTof%u%u%uSts%uextr", TofStationNo, modPos, modPosY, o))->Fill(vExtr_tof.X(), vExtr_tof.Y());
+ fHM->H2(Form("fhTof%u%u%uSts%uXYdiff", TofStationNo, modPos, modPosY, o))
+ ->Fill(closestStsHitArray[o].X, closestStsHitArray[o].Y);
+ if (!(modPos == 0 && modPosY == 0))
+ fHM->H2(Form("fhTofSts%uXYdiff_close", o))->Fill(closestStsHitArray[o].X, closestStsHitArray[o].Y);
+ }
+ }
+
+ if ((closestStsHitStations[0].Existing == true) && (closestStsHitStations[1].Existing == true)) {
+ fHM->H2(Form("fhTof%u%u%uDiffStation01", TofStationNo, modPos, modPosY))
+ ->Fill(closestStsHitStations[0].Dist, closestStsHitStations[1].Dist);
+ fHM->H2(Form("fhTof%u%u%uXDiffStation01", TofStationNo, modPos, modPosY))
+ ->Fill(closestStsHitStations[0].X, closestStsHitStations[1].X);
+ fHM->H2(Form("fhTof%u%u%uYDiffStation01", TofStationNo, modPos, modPosY))
+ ->Fill(closestStsHitStations[0].Y, closestStsHitStations[1].Y);
+
+ if (TofStationNo == 0 && modPos == 1 && modPosY == 0) {
+ Double_t Y_new = closestStsHitStations[1].Y - closestStsHitStations[0].Y;
+ fHM->H2("fhTof010YDiffStation01_pro")->Fill(closestStsHitStations[0].Y, Y_new);
+ Double_t Y_ratio = closestStsHitStations[1].Y / closestStsHitStations[0].Y;
+ fHM->H2("fhTof010YDiffStation01_ratio")->Fill(closestStsHitStations[0].Y, Y_ratio);
+ }
+ }
+
+ } // TofHit Loop
+
+
+ // Tof Tracks
+ ClosestCbmHit closestTrackStsHitStations[2];
+
+ auto nofTofTracksEv = ev->GetNofData(ECbmDataType::kTofTrack);
+ fHM->H1("fhTofTrackMulti")->Fill(nofTofTracksEv);
+
+ Int_t Multipl_2hit_Tracks = 0;
+
+ for (int j = 0; j < nofTofTracksEv; j++) {
+ auto iTofTrack = ev->GetIndex(ECbmDataType::kTofTrack, j);
+ CbmTofTracklet* tofTrack = static_cast<CbmTofTracklet*>(fTofTrack->At(iTofTrack));
+
+ if (tofTrack->GetNofHits() == 2) Multipl_2hit_Tracks++;
+
+ //Sts Loop
+ auto nofStsHitsEv = ev->GetNofData(ECbmDataType::kStsHit);
+ for (int jSts = 0; jSts < nofStsHitsEv; jSts++) {
+ auto iStsHit = ev->GetIndex(ECbmDataType::kStsHit, jSts);
+ CbmStsHit* stsHit = static_cast<CbmStsHit*>(fStsHit->At(iStsHit));
+
+ auto stsUnit = 0;
+ auto stsStation = (stsHit->GetAddress() >> 4) & 0xF;
+ auto stsUnitAdd = (stsHit->GetAddress() >> 8) & 0x4;
+ if (stsStation == 0 && stsUnitAdd == 0) { stsUnit = 0; }
+ else if (stsStation == 0 && stsUnitAdd != 0) {
+ stsUnit = 1;
+ }
+ else {
+ stsStation = 1;
+ stsUnit = 2;
+ }
+
+ auto Xdiff = tofTrack->GetFitX(stsHit->GetZ()) - stsHit->GetX();
+ auto Ydiff = tofTrack->GetFitY(stsHit->GetZ()) - stsHit->GetY();
+ Double_t distSts = TMath::Sqrt(Xdiff * Xdiff + Ydiff * Ydiff);
+
+ if (distSts < closestTrackStsHitStations[stsStation].Dist) {
+ closestTrackStsHitStations[stsStation].Existing = true;
+ closestTrackStsHitStations[stsStation].X = Xdiff;
+ closestTrackStsHitStations[stsStation].Y = Ydiff;
+ closestTrackStsHitStations[stsStation].Z = stsHit->GetZ();
+ closestTrackStsHitStations[stsStation].Dist = distSts;
+ closestTrackStsHitStations[stsStation].Indx = iStsHit;
+ closestTrackStsHitStations[stsStation].X_pos = stsHit->GetX();
+ closestTrackStsHitStations[stsStation].Y_pos = stsHit->GetY();
+ }
+ } // End Sts Loop
+
+ if ((closestTrackStsHitStations[0].Existing == true) && (closestTrackStsHitStations[1].Existing == true)
+ && tofTrack->GetNofHits() == 2) {
+ //std::cout<<"Station0: " << closestStsHitStations[0].Dist << " Station1: "<< closestStsHitStations[1].Dist << std::endl;
+ fHM->H2("fhTofTrackDiffStation01")
+ ->Fill(closestTrackStsHitStations[0].Dist, closestTrackStsHitStations[1].Dist);
+ fHM->H2("fhTofTrackXDiffStation01")->Fill(closestTrackStsHitStations[0].X, closestTrackStsHitStations[1].X);
+ fHM->H2("fhTofTrackYDiffStation01")->Fill(closestTrackStsHitStations[0].Y, closestTrackStsHitStations[1].Y);
+
+ fHM->H2("fhTofTrackXYDiffStation0")->Fill(closestTrackStsHitStations[0].X, closestTrackStsHitStations[0].Y);
+ fHM->H2("fhTofTrackXYDiffStation1")->Fill(closestTrackStsHitStations[1].X, closestTrackStsHitStations[1].Y);
+
+
+ fHM->H2("fhTofTrackXY520cm_2Hits")->Fill(tofTrack->GetFitX(520.0), tofTrack->GetFitY(520.0));
+ fHM->H2("fhTofTrackXYVertex_2Hits")->Fill(tofTrack->GetFitX(0.0), tofTrack->GetFitY(0.0));
+
+ // Time Offset Calculation
+ for (auto cntStation = 0; cntStation < 2; ++cntStation) {
+ for (auto kTrackHits = 0; kTrackHits < tofTrack->GetNofHits(); ++kTrackHits) {
+ auto trackHitInd = tofTrack->GetHitIndex(kTrackHits);
+ CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHit->At(trackHitInd));
+
+ if (tofHit == nullptr) continue;
+
+ auto diffX = tofHit->GetX() - closestTrackStsHitStations[cntStation].X_pos;
+ auto diffY = tofHit->GetY() - closestTrackStsHitStations[cntStation].Y_pos;
+ auto diffZ = tofHit->GetZ() - closestTrackStsHitStations[cntStation].Z;
+
+ auto distance = TMath::Sqrt(diffX * diffX + diffY * diffY + diffZ * diffZ);
+ Double_t TimeOfFlight = distance * tofTrack->GetTt(); // ns
+ if (tofTrack->GetTt() < 0) continue;
+ fHM->H1("fhTimeOfFlight")->Fill(TimeOfFlight);
+
+ auto TofStationNo = (tofHit->GetAddress() >> 4) & 0xf;
+ fHM->H1(Form("fhTimeOfFlight_%u", TofStationNo))->Fill(TimeOfFlight);
+
+ fHM->H1(Form("fhTimeOfFlight_%u_%u", TofStationNo, cntStation))->Fill(TimeOfFlight);
+
+ fHM->H2(Form("fhTimeOfFlight_%u_diffZ", TofStationNo))->Fill(diffZ, TimeOfFlight);
+
+ // tof hit - TimeOfFlight - stshittime
+ CbmStsHit* stsHit = static_cast<CbmStsHit*>(fStsHit->At(closestTrackStsHitStations[cntStation].Indx));
+ if (nullptr == stsHit) continue;
+
+ auto timeoffset = tofHit->GetTime() - TimeOfFlight - stsHit->GetTime();
+ fHM->H1(Form("fhTimeOffset"))->Fill(timeoffset);
+ }
+ }
+ }
+
+ fHM->H1("fhTofHitsPerTrack")->Fill(tofTrack->GetNofHits());
+ //std::cout<< tofTrack->GetT0() << " " << tofTrack->GetTt() << " " << tofTrack->GetT0()-tofTrack->GetTt() <<std::endl;
+ fHM->H1("fhTofTracksT0")->Fill(tofTrack->GetT0() - eventStartTime);
+ fHM->H1("fhTofTracksTt")->Fill(tofTrack->GetTt()); // invers velocity
+ fHM->H1("fhTofTracksVelocity")->Fill(1. / tofTrack->GetTt()); // invers velocity
+ fHM->H1("fhTofTracksTdiff")->Fill(tofTrack->GetTime() - eventStartTime);
+ fHM->H1("fhTofTracksTdiff2")->Fill(tofTrack->GetTime() - tofTrack->GetT0());
+ fHM->H1("fhTofTracksTime")->Fill(tofTrack->GetTime());
+ fHM->H1("fhTofTracksDistance")->Fill(tofTrack->GetDistance());
+ fHM->H2("fhTofTrackXY_basic")->Fill(tofTrack->GetTrackX(), tofTrack->GetTrackY());
+ Double_t inVel = 1e7 / (tofTrack->GetTt()); // in m/s
+ auto beta = inVel / TMath::C();
+ fHM->H1("fhTofBeta")->Fill(beta);
+
+ fHM->H2("fhTofTrackXYVertex")->Fill(tofTrack->GetFitX(0.0), tofTrack->GetFitY(0.0));
+ fHM->H2("fhTofTrackXY520cm")->Fill(tofTrack->GetFitX(520.0), tofTrack->GetFitY(520.0));
+ fHM->H2("fhTofTrackXY")->Fill(tofTrack->GetFitX(260.), tofTrack->GetFitY(260.));
+
+ auto nofTrackHits = tofTrack->GetNofHits();
+ for (auto kTrackHits = 0; kTrackHits < nofTrackHits; ++kTrackHits) {
+ auto trackHitInd = tofTrack->GetHitIndex(kTrackHits);
+ CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHit->At(trackHitInd));
+
+ fHM->H1("fhTofTrackHitsTimeDiff")->Fill(tofHit->GetTime() - tofTrack->GetT0());
+ fHM->H1("fhTofTrackHitsTime")->Fill(tofHit->GetTime() - eventStartTime);
+ }
+
+ Double_t closestStsR = 100000.0;
+ Int_t closestStsIndx = -1;
+
+ for (int jSts = 0; jSts < nofStsHitsEv; jSts++) {
+ auto iStsHit = ev->GetIndex(ECbmDataType::kStsHit, jSts);
+ CbmStsHit* stsHit = static_cast<CbmStsHit*>(fStsHit->At(iStsHit));
+
+ auto Xdiff = tofTrack->GetFitX(stsHit->GetZ()) - stsHit->GetX();
+ auto Ydiff = tofTrack->GetFitY(stsHit->GetZ()) - stsHit->GetY();
+
+ fHM->H2("fhTofTrackStsXYdiff")->Fill(Xdiff, Ydiff);
+
+ auto rStsTofTrack = TMath::Sqrt(Xdiff * Xdiff + Ydiff * Ydiff);
+ if (rStsTofTrack < closestStsR) {
+ closestStsR = rStsTofTrack;
+ closestStsIndx = iStsHit;
+ }
+ } // Sts Hit Loop
+
+ if (closestStsIndx > -1 && closestStsR < CutCloseDistStsTrack) {
+ CbmStsHit* closeStsHit = static_cast<CbmStsHit*>(fStsHit->At(closestStsIndx));
+ fHM->H1("fhTofTrackStsR")->Fill(closestStsR);
+
+ // 1.) calculate length from T0 to Sts Hit
+ auto stsT0Length = std::sqrt(std::pow((closeStsHit->GetX() - tofTrack->GetTrackX()), 2)
+ + std::pow((closeStsHit->GetY() - tofTrack->GetTrackY()), 2)
+ + std::pow((closeStsHit->GetZ() - 0.0), 2));
+ fHM->H1("fhTofTracksLength")->Fill(stsT0Length);
+
+ // 2.) calculate time in Tof system with length and velocity
+ auto stsTrackTime = stsT0Length * tofTrack->GetTt(); // time (of flight) in ns
+ fHM->H1("fhStsTrackTime")->Fill(stsTrackTime);
+
+ // 3.) compare to Sts hit Time
+ auto diffStsTimes = (tofTrack->GetT0() + stsTrackTime) - closeStsHit->GetTime();
+ fHM->H1("fhStsHitStsTrackTimeDiff")->Fill(diffStsTimes);
+ }
+ } // Track Loop
+
+ fHM->H1("fhMultiplicityTof2HitsTracks")->Fill(Multipl_2hit_Tracks);
+ } // Event Loop
+ } // TS loop
+
+ outFileLin->Close();
+ /// Restore old global file and folder pointer to avoid messing with FairRoot
+ gFile = oldFileLin;
+ gDirectory->cd(oldirLin->GetPath());
+
+ fRECO->Close();
+
+ std::cout << "tofHitYMin: " << tofHitYMin << " tofHitYMax: " << tofHitYMax << std::endl;
+
+ {
+ ofstream XY_differences;
+ XY_differences.open(Form("%s/XY_differences.dat", outdir));
+ // Alignments:
+ for (auto l = 0; l < 3; ++l) {
+ for (auto i = 0; i < 2; ++i) {
+ for (auto m = 0; m < 5; ++m) {
+ auto j = 0;
+ auto k = 0;
+
+ if (m == 0) {
+ j = 0;
+ k = 0;
+ }
+ if (m == 1) {
+ j = 1;
+ k = 0;
+ }
+ if (m == 2) {
+ j = 0;
+ k = 1;
+ }
+ if (m == 3) {
+ j = 1;
+ k = 1;
+ }
+ if (m == 4) {
+ j = 0;
+ k = 2;
+ }
+
+ auto x = fHM->H2(Form("fhTof%u%u%uSts%uXYdiff", i, j, k, l))->GetMean(1);
+ auto y = fHM->H2(Form("fhTof%u%u%uSts%uXYdiff", i, j, k, l))->GetMean(2);
+ std::cout << " Tof:" << i << j << k << " StS:" << l << " X:" << x << " Y:" << y << std::endl;
+ XY_differences << " Tof:" << i << j << k << " StS:" << l << " X:" << x << " Y:" << y << std::endl;
+ }
+ }
+ std::cout << std::endl;
+ XY_differences << std::endl;
+ }
+ XY_differences.close();
+ }
+
+
+ std::vector<TH1*> fitHistsX, fitHistsY;
+ std::vector<std::string> XFitstring, YFitstring;
+
+ for (auto p = 0; p < 3; ++p) {
+ ProjectXY(Form("fhTofSts%uXYdiff_close", p), fHM, &fitHistsX, &XFitstring, &fitHistsY, &YFitstring);
+ }
+
+ {
+ std::vector<std::string> strTofStsHistX, strTofStsHistY;
+ ProjectXY("fhTof010Sts0XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof001Sts0XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof011Sts0XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof002Sts0XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+
+ ProjectXY("fhTof110Sts0XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof101Sts0XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof111Sts0XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof102Sts0XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+
+
+ ProjectXY("fhTof010Sts1XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof001Sts1XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof011Sts1XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof002Sts1XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+
+ ProjectXY("fhTof110Sts1XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof101Sts1XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof111Sts1XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof102Sts1XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+
+
+ ProjectXY("fhTof010Sts2XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof001Sts2XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof011Sts2XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof002Sts2XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+
+ ProjectXY("fhTof110Sts2XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof101Sts2XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof111Sts2XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+ ProjectXY("fhTof102Sts2XYdiff", fHM, &fitHistsX, &strTofStsHistX, &fitHistsY, &strTofStsHistY);
+
+ std::cout << "--------------------------------------------------------------------" << std::endl;
+ std::cout << "fhTof*10Sts0XYdiff : " << strTofStsHistX.at(0) << " | " << strTofStsHistY.at(0) << " || "
+ << strTofStsHistX.at(4) << " | " << strTofStsHistY.at(4) << std::endl;
+ std::cout << "fhTof*01Sts0XYdiff : " << strTofStsHistX.at(1) << " | " << strTofStsHistY.at(1) << " || "
+ << strTofStsHistX.at(5) << " | " << strTofStsHistY.at(5) << std::endl;
+ std::cout << "fhTof*11Sts0XYdiff : " << strTofStsHistX.at(2) << " | " << strTofStsHistY.at(2) << " || "
+ << strTofStsHistX.at(6) << " | " << strTofStsHistY.at(6) << std::endl;
+ std::cout << "fhTof*02Sts0XYdiff : " << strTofStsHistX.at(3) << " | " << strTofStsHistY.at(3) << " || "
+ << strTofStsHistX.at(7) << " | " << strTofStsHistY.at(7) << std::endl;
+ std::cout << "--------------------------------------------------------------------" << std::endl;
+
+ std::cout << "--------------------------------------------------------------------" << std::endl;
+ std::cout << "fhTof*10Sts1XYdiff : " << strTofStsHistX.at(8) << " | " << strTofStsHistY.at(8) << " || "
+ << strTofStsHistX.at(12) << " | " << strTofStsHistY.at(12) << std::endl;
+ std::cout << "fhTof*01Sts1XYdiff : " << strTofStsHistX.at(9) << " | " << strTofStsHistY.at(9) << " || "
+ << strTofStsHistX.at(13) << " | " << strTofStsHistY.at(13) << std::endl;
+ std::cout << "fhTof*11Sts1XYdiff : " << strTofStsHistX.at(10) << " | " << strTofStsHistY.at(10) << " || "
+ << strTofStsHistX.at(14) << " | " << strTofStsHistY.at(14) << std::endl;
+ std::cout << "fhTof*02Sts1XYdiff : " << strTofStsHistX.at(11) << " | " << strTofStsHistY.at(11) << " || "
+ << strTofStsHistX.at(15) << " | " << strTofStsHistY.at(15) << std::endl;
+ std::cout << "--------------------------------------------------------------------" << std::endl;
+
+ std::cout << "--------------------------------------------------------------------" << std::endl;
+ std::cout << "fhTof*10Sts2XYdiff : " << strTofStsHistX.at(16) << " | " << strTofStsHistY.at(16) << " || "
+ << strTofStsHistX.at(20) << " | " << strTofStsHistY.at(20) << std::endl;
+ std::cout << "fhTof*01Sts2XYdiff : " << strTofStsHistX.at(17) << " | " << strTofStsHistY.at(17) << " || "
+ << strTofStsHistX.at(21) << " | " << strTofStsHistY.at(21) << std::endl;
+ std::cout << "fhTof*11Sts2XYdiff : " << strTofStsHistX.at(18) << " | " << strTofStsHistY.at(18) << " || "
+ << strTofStsHistX.at(22) << " | " << strTofStsHistY.at(22) << std::endl;
+ std::cout << "fhTof*02Sts2XYdiff : " << strTofStsHistX.at(19) << " | " << strTofStsHistY.at(19) << " || "
+ << strTofStsHistX.at(23) << " | " << strTofStsHistY.at(23) << std::endl;
+ std::cout << "--------------------------------------------------------------------" << std::endl;
+ }
+ drawHists(fHM);
+
+ std::cout << "TimeOffset Fit: " << FitTimeOffset(fHM->H1("fhTimeOffset")) << std::endl;
+
+ if (true) {
+ /// Save old global file and folder pointer to avoid messing with FairRoot
+ TFile* oldFile = gFile;
+ TDirectory* oldir = gDirectory;
+
+ std::string s = Form("%s/RecoInspectAlignHists.root", outdir);
+ TFile* outFile = new TFile(s.c_str(), "RECREATE");
+ if (outFile->IsOpen()) {
+ fHM->WriteToFile();
+ for (auto xhists = 0; xhists < fitHistsX.size(); xhists++)
+ fitHistsX.at(xhists)->Write();
+ for (auto yhists = 0; yhists < fitHistsY.size(); yhists++)
+ fitHistsY.at(yhists)->Write();
+ std::cout << "Written to Root-file \"" << s << "\" ...";
+ outFile->Close();
+ std::cout << "Done!" << std::endl;
+ }
+ /// Restore old global file and folder pointer to avoid messing with FairRoot
+ gFile = oldFile;
+ gDirectory->cd(oldir->GetPath());
+ }
+
+ std::cout << " --------------------------------------------" << std::endl;
+ std::cout << "0 | " << XFitstring.at(0) << " | " << YFitstring.at(0) << std::endl;
+ std::cout << "1 | " << XFitstring.at(1) << " | " << YFitstring.at(1) << std::endl;
+ std::cout << "2 | " << XFitstring.at(2) << " | " << YFitstring.at(2) << std::endl;
+}
+
+
+void initHists(CbmHistManager* fHM, Int_t nofEvents)
+{
+ fHM->Create1<TH1D>("fhNofTimeslices", "fhNofTimeslices;Entries", 1, 0.5, 1.5);
+
+ fHM->Create2<TH2D>("fhTofHitsXY", "fhTofHitsXY;X [cm];Y [cm];Entries", 200, -100., 100., 200, -100., 100.);
+ fHM->Create1<TH1D>("fhTofTrackMulti", "fhTofTrackMulti;Tracks/Event;Entries", 20, -0.5, 19.5);
+ fHM->Create1<TH1D>("fhTofHitsZ", "fhTofHitsZ;Z[cm];Entries", 100, 250, 300);
+ fHM->Create1<TH1D>("fhTofHitsPerTrack", "fhTofHitsPerTrack;Hits/Track;Entries", 11, -0.5, 10.5);
+ fHM->Create1<TH1D>("fhTofTracksT0", "fhTofTracksT0;T0 time [ns];Entries", 100, -50, 50);
+ fHM->Create1<TH1D>("fhTofTracksTime", "fhTofTracksTime;Track time [ns];Entries", 30000, -50, 3.0e8);
+ fHM->Create1<TH1D>("fhTofTracksTt", "fhTofTracksTt;inv. velocity [ns/cm];Entries", 100, -50, 50);
+ fHM->Create1<TH1D>("fhTofTracksVelocity", "fhTofTracksVelocity;velocity [cm/ns];Entries", 100, -50, 50);
+ fHM->Create1<TH1D>("fhTofTracksTdiff", "fhTofTracksTdiff;Time [ns];Entries", 100, -50, 50);
+ fHM->Create1<TH1D>("fhTofTracksTdiff2", "fhTofTracksTdiff2;Time [ns];Entries", 100, -50, 50);
+ fHM->Create1<TH1D>("fhTofTracksLength", "fhTofTracksLength;length [cm];Entries", 400, -50, 350);
+
+ fHM->Create1<TH1D>("fhTofTrackHitsTimeDiff", "fhTofTrackHitsTimeDiff;time [ns];Entries", 100, -50.0, 50.0);
+ fHM->Create1<TH1D>("fhTofHitsTime", "fhTofHitsTime;time [ns];Entries", 100, -50.0, 50.0);
+ fHM->Create1<TH1D>("fhTofHitsTimeLayer0", "fhTofHitsTimeLayer0;time [ns];Entries", 100, -50.0, 50.0);
+ fHM->Create1<TH1D>("fhTofHitsTimeLayer1", "fhTofHitsTimeLayer1;time [ns];Entries", 100, -50.0, 50.0);
+ fHM->Create1<TH1D>("fhTofHitsTime_001", "fhTofHitsTime_001;time [ns];Entries", 100, -50.0, 50.0);
+
+ fHM->Create1<TH1D>("fhTofTrackHitsTime", "fhTofTrackHitsTime;time [ns];Entries", 100, -50.0, 50.0);
+
+ fHM->Create1<TH1D>("fhTofBeta", "fhTofBeta;#beta;Entries", 2200, -1.1, 1.1);
+ fHM->Create1<TH1D>("fhTofTracksDistance", "fhTofTracksDistance;distance [cm];Entries", 100, -5, 5);
+
+ fHM->Create2<TH2D>("fhTofTrackXYVertex", "fhTofTrackXYVertex;X [cm];Y [cm];Entries", 200, -100, 100, 200, -100, 100);
+ fHM->Create2<TH2D>("fhTofTrackXYVertex_2Hits", "fhTofTrackXYVertex_2Hits;X [cm];Y [cm];Entries", 200, -100, 100, 200,
+ -100, 100);
+ fHM->Create2<TH2D>("fhTofTrackXY520cm", "fhTofTrackXY520cm;X [cm];Y [cm];Entries", 200, -100, 100, 200, -100, 100);
+ fHM->Create2<TH2D>("fhTofTrackXY520cm_2Hits", "fhTofTrackXY520cm_2Hits;X [cm];Y [cm];Entries", 200, -100, 100, 200,
+ -100, 100);
+ fHM->Create2<TH2D>("fhTofTrackXY", "fhTofTrackXY;X [cm];Y [cm];Entries", 200, -100, 100, 200, -100, 100);
+ fHM->Create2<TH2D>("fhTofTrackXY_basic", "fhTofTrackXY_basic;X [cm];Y [cm];Entries", 200, -100, 100, 200, -100, 100);
+
+
+ fHM->Create1<TH1D>("fhTofTrackStsR", "fhTofTrackStsR;distance [cm];Entries", 600, 0.0, 60.0);
+ fHM->Create1<TH1D>("fhStsHitStsTrackTimeDiff", "fhStsHitStsTrackTimeDiff;time diff [ns];Entries", 1200, -60.0, 60.0);
+ fHM->Create1<TH1D>("fhStsTrackTime", "fhStsTrackTime;time from track [ns];Entries", 1200, -60.0, 60.0);
+ fHM->Create1<TH1D>("fhStsHitsZ", "fhStsHitsZ; Z [cm];Entries", 400, 20.0, 60.0);
+
+ fHM->Create2<TH2D>("fhTofTrackStsXYdiff", "fhTofTrackStsXYdiff;X_tof-X_Sts [cm];Y_tof-Y_Sts [cm];Entries", 600, -2.95,
+ 2.95, 600, -2.95, 2.95);
+
+ for (auto i = 0; i < 2; ++i) {
+ for (auto j = 0; j < 2; ++j) {
+ for (auto k = 0; k < 3; ++k) {
+ if (k == 2 && j == 1) continue;
+ fHM->Create2<TH2D>(Form("fhTof%u%u%uDiffStation01", i, j, k),
+ Form("fhTof%u%u%uDiffStation01;Station 0 [cm];Station1 [cm];Entries", i, j, k), 200, -0.5,
+ 19.5, 200, -0.5, 19.5);
+ fHM->Create2<TH2D>(Form("fhTof%u%u%uXDiffStation01", i, j, k),
+ Form("fhTof%u%u%uXDiffStation01;Station 0 [cm];Station1 [cm];Entries", i, j, k), 110, -5.5,
+ 5.5, 110, -5.5, 5.5);
+ fHM->Create2<TH2D>(Form("fhTof%u%u%uYDiffStation01", i, j, k),
+ Form("fhTof%u%u%uYDiffStation01;Station 0 [cm];Station1 [cm];Entries", i, j, k), 110, -5.5,
+ 5.5, 110, -5.5, 5.5);
+ for (auto l = 0; l < 3; ++l) {
+ fHM->Create2<TH2D>(Form("fhTof%u%u%uSts%uXYdiff", i, j, k, l),
+ Form("fhTof%u%u%uSts%uXYdiff;X_tof-X_Sts [cm];Y_tof-Y_Sts [cm];Entries", i, j, k, l), 300,
+ -2.95, 2.95, 300, -2.95, 2.95);
+ fHM->Create2<TH2D>(Form("fhTof%u%u%uSts%uextr", i, j, k, l),
+ Form("fhTof%u%u%uSts%uextr;X_tof [cm];Y_tof [cm];Entries", i, j, k, l), 600, -29.5, 29.5,
+ 600, -29.5, 29.5);
+ }
+ }
+ }
+ }
+
+
+ fHM->Create2<TH2D>("fhTof010YDiffStation01_pro", "fhTof010YDiffStation01_pro;Station 0 [cm];Station1 [cm];Entries",
+ 110, -5.5, 5.5, 110, -5.5, 5.5);
+ fHM->Create2<TH2D>("fhTof010YDiffStation01_ratio",
+ "fhTof010YDiffStation01_ratio;Station 0 [cm];Station1 [cm];Entries", 110, -5.5, 5.5, 110, -5.5,
+ 5.5);
+
+ fHM->Create2<TH2D>("fhTofTrackDiffStation01", "fhTofTrackDiffStation01;Station 0 [cm];Station1 [cm];Entries", 400,
+ -0.5, 39.5, 400, -0.5, 39.5);
+ fHM->Create2<TH2D>("fhTofTrackXDiffStation01", "fhTofTrackXDiffStation01;Station 0 [cm];Station1 [cm];Entries", 790,
+ -39.5, 39.5, 790, -39.5, 39.5);
+ fHM->Create2<TH2D>("fhTofTrackYDiffStation01", "fhTofTrackYDiffStation01;Station 0 [cm];Station1 [cm];Entries", 790,
+ -39.5, 39.5, 790, -39.5, 39.5);
+ fHM->Create2<TH2D>("fhTofTrackXYDiffStation0", "fhTofTrackXYDiffStation0;Xdiff [cm];Ydiff [cm];Entries", 790, -79.5,
+ 79.5, 790, -79.5, 79.5);
+ fHM->Create2<TH2D>("fhTofTrackXYDiffStation1", "fhTofTrackXYDiffStation1;Xdiff [cm];Ydiff [cm];Entries", 790, -79.5,
+ 79.5, 790, -79.5, 79.5);
+
+
+ fHM->Create2<TH2D>("fhTofSts0XYdiff", "fhTofSts0XYdiff;X_tof-X_Sts [cm];Y_tof-Y_Sts [cm];Entries", 600, -2.95, 2.95,
+ 600, -2.95, 2.95);
+ fHM->Create2<TH2D>("fhTofSts1XYdiff", "fhTofSts1XYdiff;X_tof-X_Sts [cm];Y_tof-Y_Sts [cm];Entries", 600, -2.95, 2.95,
+ 600, -2.95, 2.95);
+ fHM->Create2<TH2D>("fhTofSts2XYdiff", "fhTofSts2XYdiff;X_tof-X_Sts [cm];Y_tof-Y_Sts [cm];Entries", 600, -2.95, 2.95,
+ 600, -2.95, 2.95);
+
+ fHM->Create2<TH2D>("fhTofSts0XYdiff_close", "fhTofSts0XYdiff_close;X_tof-X_Sts [cm];Y_tof-Y_Sts [cm];Entries", 600,
+ -2.95, 2.95, 600, -2.95, 2.95);
+ fHM->Create2<TH2D>("fhTofSts1XYdiff_close", "fhTofSts1XYdiff_close;X_tof-X_Sts [cm];Y_tof-Y_Sts [cm];Entries", 600,
+ -2.95, 2.95, 600, -2.95, 2.95);
+ fHM->Create2<TH2D>("fhTofSts2XYdiff_close", "fhTofSts2XYdiff_close;X_tof-X_Sts [cm];Y_tof-Y_Sts [cm];Entries", 600,
+ -2.95, 2.95, 600, -2.95, 2.95);
+
+ fHM->Create2<TH2D>("fhTofSts0XY_8.5cm", "fhTofSts0XY_8.5cm;X [cm];Y [cm];Entries", 100, -1.0, 1.0, 100, -2.5, 2.5);
+ fHM->Create2<TH2D>("fhTofSts1XY_8.5cm", "fhTofSts1XY_8.5cm;X [cm];Y [cm];Entries", 100, -1.0, 1.0, 100, -2.5, 2.5);
+ fHM->Create2<TH2D>("fhTofSts2XY_8.5cm", "fhTofSts2XY_8.5cm;X [cm];Y [cm];Entries", 100, -1.0, 1.0, 100, -2.5, 2.5);
+
+
+ fHM->Create2<TH2D>("fhStsUnit0Hits", "fhStsUnit0Hits;X_Sts [cm];Y_Sts [cm];Entries", 600, -30, 30, 600, -30, 30);
+ fHM->Create2<TH2D>("fhStsUnit1Hits", "fhStsUnit1Hits;X_Sts [cm];Y_Sts [cm];Entries", 600, -30, 30, 600, -30, 30);
+ fHM->Create2<TH2D>("fhStsUnit2Hits", "fhStsUnit2Hits;X_Sts [cm];Y_Sts [cm];Entries", 600, -30, 30, 600, -30, 30);
+
+ fHM->Create2<TH2D>("fhStsStation0Hits", "fhStsStation0Hits;X_Sts [cm];Y_Sts [cm];Entries", 600, -30, 30, 600, -30,
+ 30);
+
+ fHM->Create2<TH2D>("fhTofHitsXY_001", "fhTofHitsXY_001;X_tof [cm];Y_tof [cm];Entries", 80, -20, 20, 80, -20, 20);
+
+ fHM->Create1<TH1D>("fhLinearFitZX", "fhLinearFitZX;X@0 [cm];Entries", 400, -2, 2);
+ fHM->Create1<TH1D>("fhLinearFitZY", "fhLinearFitZY;Y@0 [cm];Entries", 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY", "fhLinearFitXY;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create1<TH1D>("fhLinearFitChi2X", "fhLinearFitChi2X;Chi2 X-Fit;Entries", 200, 0, 2);
+ fHM->Create1<TH1D>("fhLinearFitChi2Y", "fhLinearFitChi2Y;Chi2 Y-Fit;Entries", 200, 0, 2);
+
+ fHM->Create2<TH2D>("fhLinearFitXY_.25", "fhLinearFitXY_.25;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_.50", "fhLinearFitXY_.50;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_.75", "fhLinearFitXY_.75;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_1", "fhLinearFitXY_2;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_2", "fhLinearFitXY_2;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_4", "fhLinearFitXY_4;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_6", "fhLinearFitXY_6;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+
+ fHM->Create2<TH2D>("fhLinearFitXY_-.5", "fhLinearFitXY_-.5;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_-1", "fhLinearFitXY_-1;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_-2", "fhLinearFitXY_-2;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_-4", "fhLinearFitXY_-4;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+ fHM->Create2<TH2D>("fhLinearFitXY_-6", "fhLinearFitXY_-6;X@0 [cm];Y@0 [cm];Entries", 400, -2, 2, 400, -2, 2);
+
+ fHM->Create1<TH1D>("fhTimeOffset", "fhTimeOffset;time [ns];Entries", 800, -40.0, 40.0);
+ fHM->Create1<TH1D>("fhTimeOfFlight", "fhTimeOfFlight;time [ns];Entries", 800, -40.0, 40.0);
+ fHM->Create1<TH1D>("fhTimeOfFlight_0", "fhTimeOfFlight_0;time [ns];Entries", 800, -40.0, 40.0);
+ fHM->Create1<TH1D>("fhTimeOfFlight_1", "fhTimeOfFlight_1;time [ns];Entries", 800, -40.0, 40.0);
+
+ fHM->Create1<TH1D>("fhTimeOfFlight_0_0", "fhTimeOfFlight_0_0;time [ns];Entries", 800, -40.0, 40.0);
+ fHM->Create1<TH1D>("fhTimeOfFlight_0_1", "fhTimeOfFlight_0_1;time [ns];Entries", 800, -40.0, 40.0);
+ fHM->Create1<TH1D>("fhTimeOfFlight_1_0", "fhTimeOfFlight_1_0;time [ns];Entries", 800, -40.0, 40.0);
+ fHM->Create1<TH1D>("fhTimeOfFlight_1_1", "fhTimeOfFlight_1_1;time [ns];Entries", 800, -40.0, 40.0);
+
+ fHM->Create2<TH2D>("fhTimeOfFlight_0_diffZ", "fhTimeOfFlight_0_diffZ;time [ns];Entries", 300, 0., 300., 800, -40.0,
+ 40.0);
+ fHM->Create2<TH2D>("fhTimeOfFlight_1_diffZ", "fhTimeOfFlight_1_diffZ;time [ns];Entries", 300, 0., 300., 800, -40.0,
+ 40.0);
+
+ fHM->Create1<TH1D>("fhSimpleTrackTofTime", "fhSimpleTrackTofTime;time [ns];Entries", 400, -2.0, 2.0);
+ fHM->Create1<TH1D>("fhSimpleTrackTofDist", "fhSimpleTrackTofDist;distance [cm];Entries", 150, 10.0, 25.0);
+ fHM->Create1<TH1D>("fhSimpleTrackTofBeta", "fhSimpleTrackTofBeta;#beta;Entries", 600, -3.0, 3.0);
+
+ fHM->Create1<TH1D>("fhMultiplicityStsHits", "fhMultiplicityStsHits;Hits/Event;Entries", 20, 0.0, 20.0);
+ fHM->Create1<TH1D>("fhMultiplicityTofHits", "fhMultiplicityTofHits;Hits/Event;Entries", 20, 0.0, 20.0);
+ fHM->Create1<TH1D>("fhMultiplicityTof2HitsTracks", "fhMultiplicityTof2HitsTracks;2-hit Tracks/Event;Entries", 20, 0.0,
+ 20.0);
+
+ auto EventsPerTSBinning = (nofEvents < 200) ? nofEvents : 200;
+ fHM->Create1<TH1D>("fhEventsPerTS", "fhEventsPerTS;TS;#Events", EventsPerTSBinning, 0, EventsPerTSBinning);
+}
+
+void drawHists(CbmHistManager* fHM)
+{
+
+ {
+ TCanvas* c = fHM->CreateCanvas("EventsPerTS", "EventsPerTS", 2000, 2000);
+ fHM->H1("fhEventsPerTS")->Draw("HIST");
+
+ c->SaveAs(Form("%s/EventsPerTS.png", outdir));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("StsUnit0_TofXYDiff", "StsUnit0_TofXYDiff", 2000, 2000);
+ c->Divide(4, 5);
+ c->cd(1);
+ fHM->H2("fhTof000Sts0XYdiff")->Draw("colz");
+ c->cd(3);
+ fHM->H2("fhTof100Sts0XYdiff")->Draw("colz");
+
+ c->cd(6);
+ fHM->H2("fhTof010Sts0XYdiff")->Draw("colz");
+ c->cd(8);
+ fHM->H2("fhTof110Sts0XYdiff")->Draw("colz");
+
+ c->cd(9);
+ fHM->H2("fhTof001Sts0XYdiff")->Draw("colz");
+ c->cd(11);
+ fHM->H2("fhTof101Sts0XYdiff")->Draw("colz");
+
+ c->cd(14);
+ fHM->H2("fhTof011Sts0XYdiff")->Draw("colz");
+ c->cd(16);
+ fHM->H2("fhTof111Sts0XYdiff")->Draw("colz");
+
+ c->cd(17);
+ fHM->H2("fhTof002Sts0XYdiff")->Draw("colz");
+
+ c->cd(19);
+ fHM->H2("fhTof102Sts0XYdiff")->Draw("colz");
+
+ c->SaveAs(Form("%s/StsUnit0_TofXYDiff.png", outdir));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("StsUnit1_TofXYDiff", "StsUnit1_TofXYDiff", 2000, 2000);
+ c->Divide(4, 5);
+ c->cd(1);
+ fHM->H2("fhTof000Sts1XYdiff")->Draw("colz");
+ c->cd(3);
+ fHM->H2("fhTof100Sts1XYdiff")->Draw("colz");
+
+ c->cd(6);
+ fHM->H2("fhTof010Sts1XYdiff")->Draw("colz");
+ c->cd(8);
+ fHM->H2("fhTof110Sts1XYdiff")->Draw("colz");
+
+ c->cd(9);
+ fHM->H2("fhTof001Sts1XYdiff")->Draw("colz");
+ c->cd(11);
+ fHM->H2("fhTof101Sts1XYdiff")->Draw("colz");
+
+ c->cd(14);
+ fHM->H2("fhTof011Sts1XYdiff")->Draw("colz");
+ c->cd(16);
+ fHM->H2("fhTof111Sts1XYdiff")->Draw("colz");
+
+ c->cd(17);
+ fHM->H2("fhTof002Sts1XYdiff")->Draw("colz");
+
+ c->cd(19);
+ fHM->H2("fhTof102Sts1XYdiff")->Draw("colz");
+
+ c->SaveAs(Form("%s/StsUnit1_TofXYDiff.png", outdir));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("StsUnit2_TofXYDiff", "StsUnit2_TofXYDiff", 2000, 2000);
+ c->Divide(4, 5);
+ c->cd(1);
+ fHM->H2("fhTof000Sts2XYdiff")->Draw("colz");
+ c->cd(3);
+ fHM->H2("fhTof100Sts2XYdiff")->Draw("colz");
+
+ c->cd(6);
+ fHM->H2("fhTof010Sts2XYdiff")->Draw("colz");
+ c->cd(8);
+ fHM->H2("fhTof110Sts2XYdiff")->Draw("colz");
+
+ c->cd(9);
+ fHM->H2("fhTof001Sts2XYdiff")->Draw("colz");
+ c->cd(11);
+ fHM->H2("fhTof101Sts2XYdiff")->Draw("colz");
+
+ c->cd(14);
+ fHM->H2("fhTof011Sts2XYdiff")->Draw("colz");
+ c->cd(16);
+ fHM->H2("fhTof111Sts2XYdiff")->Draw("colz");
+
+ c->cd(17);
+ fHM->H2("fhTof002Sts2XYdiff")->Draw("colz");
+
+ c->cd(19);
+ fHM->H2("fhTof102Sts2XYdiff")->Draw("colz");
+
+ c->SaveAs(Form("%s/StsUnit2_TofXYDiff.png", outdir));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("StsUnits_TofXYDiff", "StsUnits_TofXYDiff", 3000, 1000);
+ c->Divide(3, 1);
+ c->cd(1);
+ fHM->H2("fhTofSts0XYdiff")->Draw("colz");
+ c->cd(2);
+ fHM->H2("fhTofSts1XYdiff")->Draw("colz");
+ c->cd(3);
+ fHM->H2("fhTofSts2XYdiff")->Draw("colz");
+
+ c->SaveAs(Form("%s/StsUnits_TofXYDiff.png", outdir));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("Multiplicities", "Multiplicities", 3000, 1000);
+ c->Divide(3, 1);
+ c->cd(1);
+ fHM->H1("fhMultiplicityStsHits")->Draw("Hist");
+ c->cd(2);
+ fHM->H1("fhMultiplicityTofHits")->Draw("Hist");
+ c->cd(3);
+ fHM->H1("fhMultiplicityTof2HitsTracks")->Draw("Hist");
+
+ c->SaveAs(Form("%s/Multiplicities_Sts_Tof.png", outdir));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("StsUnits_TofXYDiff_Close", "StsUnits_TofXYDiff_Close", 3000, 1000);
+ c->Divide(3, 1);
+ c->cd(1);
+ fHM->H2("fhTofSts0XYdiff_close")->Draw("colz");
+ c->cd(2);
+ fHM->H2("fhTofSts1XYdiff_close")->Draw("colz");
+ c->cd(3);
+ fHM->H2("fhTofSts2XYdiff_close")->Draw("colz");
+
+ c->SaveAs(Form("%s/StsUnits_TofXYDiff_Close.png", outdir));
+ }
+
+ {
+ TCanvas* c = fHM->CreateCanvas("StsZ", "StsZ", 1000, 1000);
+ fHM->H1("fhStsHitsZ")->Draw("HIST");
+
+ c->SaveAs(Form("%s/StsHitsZ.png", outdir));
+ }
+}
+
+TVector3 extrpolate(CbmTofHit* tofHit, TVector3* vertex, Double_t Z)
+{
+ TVector3 extVec(0, 0, 0);
+
+ Double_t factor = (Z - vertex->Z()) / (tofHit->GetZ() - vertex->Z());
+ Double_t x = vertex->X() + factor * (tofHit->GetX() - vertex->X());
+ Double_t y = vertex->Y() + factor * (tofHit->GetY() - vertex->Y());
+ extVec.SetXYZ(x, y, Z);
+
+ return extVec;
+}
+
+void ProjectXY(std::string name, CbmHistManager* fHM, std::vector<TH1*>* HistoX, std::vector<std::string>* stringX,
+ std::vector<TH1*>* HistoY, std::vector<std::string>* stringY)
+{
+ ProjectAxis(name, fHM, 'X', HistoX, stringX);
+ ProjectAxis(name, fHM, 'Y', HistoY, stringY);
+}
+
+void ProjectAxis(std::string name, CbmHistManager* fHM, char axis, std::vector<TH1*>* Histo,
+ std::vector<std::string>* string)
+{
+ std::cout << "start Fitting procedure" << std::endl;
+
+ TH1* hist = nullptr;
+ if (axis == 'x' || axis == 'X') hist = fHM->H2(name)->ProjectionX();
+ if (axis == 'y' || axis == 'Y') hist = fHM->H2(name)->ProjectionY();
+
+ // Find Peak position:
+ Int_t binmax = hist->GetMaximumBin();
+ Double_t max_x = hist->GetXaxis()->GetBinCenter(binmax);
+
+ TF1* g1s = new TF1("gS", "gaus(0)", max_x - 0.3, max_x + 0.3);
+ TF1* g1b = nullptr;
+
+ if (hist->GetMean() > max_x) { g1b = new TF1("gB", "gaus(0)", max_x + 1, 3); }
+ else {
+ g1b = new TF1("gB", "gaus(0)", -3, max_x - 1);
+ }
+ TF1* total1 = new TF1("total", "gaus(0)+gaus(3)", -3, 3);
+
+ Double_t par[6];
+ Double_t partotal[6];
+ hist->Fit(g1s, "QR");
+ hist->Fit(g1b, "QR+");
+
+ g1s->GetParameters(&par[0]);
+ g1b->GetParameters(&par[3]);
+ total1->SetParameters(par);
+
+ hist->Fit(total1, "QR+");
+ total1->GetParameters(&partotal[0]);
+
+ // std::cout << " Mean1: " << partotal[1] << " DEV: " << partotal[2] << " Mean2: " << partotal[4]
+ // << " DEV: " << partotal[5] << std::endl;
+ // std::cout << " MeanSig: " << par[1] << " DEV: " << par[2] << " MeanBG: " << par[4]
+ // << " DEV: " << par[5] << std::endl;
+ if (std::abs(partotal[2]) < std::abs(partotal[5])) {
+ string->emplace_back(Form("%+.4f +- %.4f cm", partotal[1], partotal[2]));
+ }
+ else {
+ string->emplace_back(Form("%+.4f +- %.4f cm", partotal[4], partotal[5]));
+ }
+
+ TCanvas* c = nullptr;
+ if (axis == 'x' || axis == 'X')
+ c = fHM->CreateCanvas(Form("%s_projX", name.c_str()), Form("%s_projX", name.c_str()), 2000, 2000);
+ if (axis == 'y' || axis == 'Y')
+ c = fHM->CreateCanvas(Form("%s_projY", name.c_str()), Form("%s_projY", name.c_str()), 2000, 2000);
+ hist->Draw("HIST");
+ total1->Draw("SAME");
+
+ if (axis == 'x' || axis == 'X') c->SaveAs(Form("%s/%s_projX.png", outdir, name.c_str()));
+ if (axis == 'y' || axis == 'Y') c->SaveAs(Form("%s/%s_projY.png", outdir, name.c_str()));
+
+ Histo->push_back(hist);
+}
+
+
+std::string FitTimeOffset(TH1* hist)
+{
+ TF1* g1s = new TF1("gS", "gaus(0)", -20, -4);
+ TF1* g1b = new TF1("gB", "gaus(0)", -40, -24);
+
+ TF1* total1 = new TF1("total", "gaus(0)+gaus(3)", -40, 40);
+
+ Double_t par[6];
+ Double_t partotal[6];
+ hist->Fit(g1s, "QR");
+ hist->Fit(g1b, "QR+");
+
+ g1s->GetParameters(&par[0]);
+ g1b->GetParameters(&par[3]);
+ total1->SetParameters(par);
+
+ hist->Fit(total1, "QR+");
+ total1->GetParameters(&partotal[0]);
+
+ //std::string FitResult = Form("%+.4f +- %.4f ns", partotal[1], partotal[2]);
+
+ TCanvas* c = new TCanvas("TimeOffset_Fit", "TimeOffset_Fit", 2000, 2000);
+ hist->Draw("HIST");
+ total1->Draw("SAME");
+ c->SaveAs(Form("%s/TimeOffset_Fit.png", outdir));
+
+ return Form("%+.4f +- %.4f ns", partotal[1], partotal[2]);
+}
diff --git a/reco/detectors/tof/CbmTofFindTracks.cxx b/reco/detectors/tof/CbmTofFindTracks.cxx
index 2b441568b1..9ffac3e5e7 100644
--- a/reco/detectors/tof/CbmTofFindTracks.cxx
+++ b/reco/detectors/tof/CbmTofFindTracks.cxx
@@ -990,6 +990,7 @@ void CbmTofFindTracks::Exec(Option_t* opt)
for (Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) {
CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
LOG(debug) << "Process event " << iEvent << " with " << tEvent->GetNofData(ECbmDataType::kTofHit) << " hits";
+ if ((tEvent->GetNofData(ECbmDataType::kTofHit)) < 1) continue;
if (fTofHitArray) fTofHitArray->Clear("C");
Int_t iNbHits = 0;
--
GitLab