CbmRoot
Loading...
Searching...
No Matches
CbmTaskTofClusterizerParWrite.cxx
Go to the documentation of this file.
1/* Copyright (C) 2022-2025 Facility for Antiproton and Ion Research in Europe, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Dominik Smith [committer], Pierre-Alain Loizeau, Norbert Herrmann */
4
5#include "CbmTaskTofClusterizerParWrite.h"
6
7// TOF Classes and includes
8//#include "CbmRunDatabaseContainer.h"
9#include "CbmTofAddress.h" // in cbmdata/tof
10#include "CbmTofCell.h" // in tof/TofData
11#include "CbmTofCreateDigiPar.h" // in tof/TofTools
12#include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
13#include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
14#include "CbmTofDetectorId_v21a.h" // in cbmdata/tof
15#include "CbmTofDigiBdfPar.h" // in tof/TofParam
16#include "CbmTofDigiPar.h" // in tof/TofParam
17#include "CbmTofGeoHandler.h" // in tof/TofTools
18
19// FAIR classes and includes
20#include "CbmYaml.h"
21#include "FairRootFileSink.h"
22#include "FairRootManager.h"
23#include "FairRunAna.h"
24#include "FairRuntimeDb.h"
25#include "algo/detectors/bmon/CalibrateSetup.h"
26#include "algo/detectors/bmon/HitfindSetup.h"
27#include "algo/detectors/tof/CalibrateSetup.h"
28#include "algo/detectors/tof/HitfindSetup.h"
29
30#include <Logger.h>
31
32// ROOT Classes and includes
33#include "TClonesArray.h"
34#include "TGeoManager.h"
35#include "TGeoPhysicalNode.h"
36#include "TH2.h"
37#include "TProfile.h"
38#include "TROOT.h"
39#include "TStopwatch.h"
40#include "TVector3.h"
41
42// C++ Classes and includes
43#include <iomanip>
44#include <vector>
45
46
47/************************************************************************************/
52
53CbmTaskTofClusterizerParWrite::CbmTaskTofClusterizerParWrite(const char* name, int32_t verbose, bool /*writeDataInOut*/)
54 : FairTask(TString(name), verbose)
55 , fTofId(NULL)
56 , fDigiPar(NULL)
57 , fDigiBdfPar(NULL)
58 , fvCPTOff()
59 , fvCPTotGain()
60 , fvCPTotOff()
61 , fvCPWalk()
62 , fvCPTOffY()
63 , fvCPTOffYBinWidth()
64 , fvCPTOffYRange()
65 , fvDeadStrips()
66 , fCalMode(0)
67 , fPosYMaxScal(1.)
68 , fTotMax(100.)
69 , fTotMin(0.)
70 , fTotOff(0.)
71 , fTotMean(0.)
72 , fMaxTimeDist(1.)
73 , fdChannelDeadtime(0.)
74 , fCalParFileName("")
75 , fdTOTMax(50.)
76 , fdTOTMin(0.)
77 , fdTTotMean(2.)
78 , fdMaxTimeDist(0.)
79 , fdMaxSpaceDist(0.)
80 , fdModifySigvel(1.0)
81 , fbSwapChannelSides(false)
82{
83}
84
85CbmTaskTofClusterizerParWrite::~CbmTaskTofClusterizerParWrite() {}
86
87/************************************************************************************/
88// FairTasks inherited functions
89InitStatus CbmTaskTofClusterizerParWrite::Init()
90{
91 LOG(info) << "CbmTaskTofClusterizerParWrite initializing... expect Digis in ns units! ";
92 if (false == InitParameters()) return kFATAL;
93 if (false == InitCalibParameter()) return kFATAL;
94 if (false == InitAlgosTof()) return kFATAL;
95 if (false == InitAlgosBmon()) return kFATAL;
96 return kSUCCESS;
97}
98
99void CbmTaskTofClusterizerParWrite::SetParContainers()
100{
101 LOG(info) << "=> Get the digi parameters for tof";
102 FairRunAna* ana = FairRunAna::Instance();
103 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
104
105 fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
106
107 LOG(info) << "found " << fDigiPar->GetNrOfModules() << " cells ";
108 fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
109}
110
111
112/************************************************************************************/
113bool CbmTaskTofClusterizerParWrite::InitParameters()
114{
115
116 // Initialize the TOF GeoHandler
117 bool isSimulation = false;
118 LOG(info) << "CbmTaskTofClusterizerParWrite::InitParameters - Geometry, Mapping, ... ??";
119
120 // Get Base Container
121 FairRun* ana = FairRun::Instance();
122 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
123
124 CbmTofGeoHandler geoHandler;
125 int32_t iGeoVersion = geoHandler.Init(isSimulation);
126 if (k14a > iGeoVersion) {
127 LOG(error) << "CbmTaskTofClusterizerParWrite::InitParameters => Only compatible "
128 "with geometries after v14a !!!";
129 return false;
130 }
131 if (iGeoVersion == k14a)
132 fTofId = new CbmTofDetectorId_v14a();
133 else
134 fTofId = new CbmTofDetectorId_v21a();
135
136 // create digitization parameters from geometry file
137 CbmTofCreateDigiPar* tofDigiPar = new CbmTofCreateDigiPar("TOF Digi Producer", "TOF task");
138 LOG(info) << "Create DigiPar ";
139 tofDigiPar->Init();
140
141 fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
142 if (0 == fDigiPar) {
143 LOG(error) << "CbmTaskTofClusterizerParWrite::InitParameters => Could not obtain "
144 "the CbmTofDigiPar ";
145 return false;
146 }
147
148 fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
149 if (0 == fDigiBdfPar) {
150 LOG(error) << "CbmTaskTofClusterizerParWrite::InitParameters => Could not obtain "
151 "the CbmTofDigiBdfPar ";
152 return false;
153 }
154
155 rtdb->initContainers(ana->GetRunId());
156
157 LOG(info) << "CbmTaskTofClusterizerParWrite::InitParameter: currently " << fDigiPar->GetNrOfModules()
158 << " digi cells ";
159
160 fdMaxTimeDist = fDigiBdfPar->GetMaxTimeDist(); // in ns
161 fdMaxSpaceDist = fDigiBdfPar->GetMaxDistAlongCh(); // in cm
162
163 if (fMaxTimeDist != fdMaxTimeDist) {
164 fdMaxTimeDist = fMaxTimeDist; // modify default
165 fdMaxSpaceDist = fdMaxTimeDist * fDigiBdfPar->GetSignalSpeed()
166 * 0.5; // cut consistently on positions (with default signal velocity)
167 }
168
169 LOG(info) << " BuildCluster with MaxTimeDist " << fdMaxTimeDist << ", MaxSpaceDist " << fdMaxSpaceDist;
170
171 fvDeadStrips.resize(fDigiBdfPar->GetNbDet());
172 return true;
173}
174/************************************************************************************/
175bool CbmTaskTofClusterizerParWrite::InitCalibParameter()
176{
177 // dimension and initialize calib parameter
178 int32_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
179
180 if (fTotMean != 0.) fdTTotMean = fTotMean; // adjust target mean for TOT
181
182 fvCPTOff.resize(iNbSmTypes);
183 fvCPTotGain.resize(iNbSmTypes);
184 fvCPTotOff.resize(iNbSmTypes);
185 fvCPWalk.resize(iNbSmTypes);
186 fvCPTOffY.resize(iNbSmTypes);
187 fvCPTOffYBinWidth.resize(iNbSmTypes);
188 fvCPTOffYRange.resize(iNbSmTypes);
189 for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
190 int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
191 int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
192 fvCPTOff[iSmType].resize(iNbSm * iNbRpc);
193 fvCPTotGain[iSmType].resize(iNbSm * iNbRpc);
194 fvCPTotOff[iSmType].resize(iNbSm * iNbRpc);
195 fvCPWalk[iSmType].resize(iNbSm * iNbRpc);
196 fvCPTOffY[iSmType].resize(iNbSm * iNbRpc);
197 fvCPTOffYBinWidth[iSmType].resize(iNbSm * iNbRpc);
198 fvCPTOffYRange[iSmType].resize(iNbSm * iNbRpc);
199 for (int32_t iSm = 0; iSm < iNbSm; iSm++) {
200 for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
201
202 fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc] = 0.; // initialize
203 fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc] = 0.; // initialize
204
205 /* D.Smith 23.8.23: For testing hit time calibration. Please remove when done.
206 fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc] = 1.; // initialize
207 fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc] = 200.; // initialize
208 fvCPTOffY[iSmType][iSm * iNbRpc + iRpc] = std::vector<double>(10000, 1000.);
209 for( size_t i = 0; i < fvCPTOffY[iSmType][iSm * iNbRpc + iRpc].size(); i++ )
210 {
211 fvCPTOffY[iSmType][iSm * iNbRpc + iRpc][i] = 1000.+500.*i;
212 }
213*/
214 int32_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
215 fvCPTOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
216 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
217 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
218 fvCPWalk[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
219 int32_t nbSide = 2 - fDigiBdfPar->GetChanType(iSmType, iRpc);
220 for (int32_t iCh = 0; iCh < iNbChan; iCh++) {
221 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
222 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
223 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
224 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
225 for (int32_t iSide = 0; iSide < nbSide; iSide++) {
226 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 0.; //initialize
227 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 1.; //initialize
228 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 0.; //initialize
229 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide].resize(nbClWalkBinX);
230 for (int32_t iWx = 0; iWx < nbClWalkBinX; iWx++) {
231 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide][iWx] = 0.;
232 }
233 }
234 }
235 }
236 }
237 }
238 LOG(info) << "CbmTaskTofClusterizerParWrite::InitCalibParameter: defaults set";
239
240 if (fCalMode <= 0) {
241 return true;
242 } // Skip calibration from histograms in mode zero
243
245 // <= To prevent histos from being sucked in by the param file of the TRootManager!
246 TFile* oldFile = gFile;
247 TDirectory* oldDir = gDirectory;
248
249 LOG(info) << "CbmTaskTofClusterizerParWrite::InitCalibParameter: read histos from "
250 << "file " << fCalParFileName;
251
252 // read parameter from histos
253 if (fCalParFileName.IsNull()) return true;
254
255 TFile* calParFile = new TFile(fCalParFileName, "READ");
256 if (NULL == calParFile) {
257 if (fCalMode % 10 == 9) { //modify reference file name
258 int iCalMode = fCalParFileName.Index("tofClust") - 3;
259 fCalParFileName(iCalMode) = '3';
260 LOG(info) << "Modified CalFileName = " << fCalParFileName;
261 calParFile = new TFile(fCalParFileName, "update");
262 if (NULL == calParFile)
263 LOG(fatal) << "CbmTaskTofClusterizerParWrite::InitCalibParameter: "
264 << "file " << fCalParFileName << " does not exist!";
265
266 return true;
267 }
268 LOG(fatal) << "Calibration parameter file not existing!";
269 }
270
271 for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
272 int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
273 int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
274 TProfile* hSvel = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iSmType));
275 for (int32_t iSm = 0; iSm < iNbSm; iSm++)
276 for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
277
278 std::vector<std::vector<double>>& vCPTotGain = fvCPTotGain[iSmType][iSm * iNbRpc + iRpc];
279 std::vector<std::vector<double>>& vCPTOff = fvCPTOff[iSmType][iSm * iNbRpc + iRpc];
280 std::vector<std::vector<double>>& vCPTotOff = fvCPTotOff[iSmType][iSm * iNbRpc + iRpc];
281 std::vector<std::vector<std::vector<double>>>& vCPWalk = fvCPWalk[iSmType][iSm * iNbRpc + iRpc];
282
283 std::vector<double>& vCPTOffY = fvCPTOffY[iSmType][iSm * iNbRpc + iRpc];
284 double& vCPTOffYBinWidth = fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc];
285 double& vCPTOffYRange = fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc];
286
287 // update default parameter
288 if (NULL != hSvel) {
289 double Vscal = hSvel->GetBinContent(iSm * iNbRpc + iRpc + 1);
290 if (Vscal == 0.) Vscal = 1.;
291 Vscal *= fdModifySigvel; //1.03; // testing the effect of wrong signal velocity, FIXME
292 fDigiBdfPar->SetSigVel(iSmType, iSm, iRpc, fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * Vscal);
293 LOG(info) << "Modify " << iSmType << iSm << iRpc << " Svel by " << Vscal << " to "
294 << fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
295 }
296 TH2D* htempPos_pfx =
297 (TH2D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx", iSmType, iSm, iRpc));
298 TH2D* htempTOff_pfx =
299 (TH2D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx", iSmType, iSm, iRpc));
300 TH1D* htempTot_Mean =
301 (TH1D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Mean", iSmType, iSm, iRpc));
302 TH1D* htempTot_Off =
303 (TH1D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Off", iSmType, iSm, iRpc));
304 if (NULL != htempPos_pfx && NULL != htempTOff_pfx && NULL != htempTot_Mean && NULL != htempTot_Off) {
305 int32_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
306 //int32_t iNbinTot = htempTot_Mean->GetNbinsX();//not used any more
307 for (int32_t iCh = 0; iCh < iNbCh; iCh++) {
308 for (int32_t iSide = 0; iSide < 2; iSide++) {
309 double TotMean = htempTot_Mean->GetBinContent(iCh * 2 + 1 + iSide); //nh +1 empirical(?)
310 if (0.001 < TotMean) {
311 vCPTotGain[iCh][iSide] *= fdTTotMean / TotMean;
312 }
313 vCPTotOff[iCh][iSide] = htempTot_Off->GetBinContent(iCh * 2 + 1 + iSide);
314 }
315 double YMean = ((TProfile*) htempPos_pfx)->GetBinContent(iCh + 1); //nh +1 empirical(?)
316 double TMean = ((TProfile*) htempTOff_pfx)->GetBinContent(iCh + 1);
317 if (std::isnan(YMean) || std::isnan(TMean)) {
318 LOG(warn) << "Invalid calibration for TSRC " << iSmType << iSm << iRpc << iCh << ", use default!";
319 continue;
320 }
321 double dTYOff = YMean / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
322 if (5 == iSmType || 8 == iSmType) dTYOff = 0.; // no valid Y positon for pad counters
323 vCPTOff[iCh][0] += -dTYOff + TMean;
324 vCPTOff[iCh][1] += +dTYOff + TMean;
325
326 if (5 == iSmType || 8 == iSmType) { // for PAD counters
327 vCPTOff[iCh][1] = vCPTOff[iCh][0];
328 vCPTotGain[iCh][1] = vCPTotGain[iCh][0];
329 vCPTotOff[iCh][1] = vCPTotOff[iCh][0];
330 }
331 //if(iSmType==0 && iSm==4 && iRpc==2 && iCh==26)
332 if (iSmType == 0 && iSm == 2 && iRpc == 0)
333 //if (iSmType == 9 && iSm == 0 && iRpc == 0 && iCh == 10) // select specific channel
334 LOG(info) << "InitCalibParameter:"
335 << " TSRC " << iSmType << iSm << iRpc << iCh
336 << Form(": YMean %6.3f, TYOff %6.3f, TMean %6.3f", YMean, dTYOff, TMean) << " -> "
337 << Form(" TOff %f, %f, TotG %f, %f ", vCPTOff[iCh][0], vCPTOff[iCh][1], vCPTotGain[iCh][0],
338 vCPTotGain[iCh][1]);
339
340 TH1D* htempWalk0 = (TH1D*) gDirectory->FindObjectAny(
341 Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh));
342 TH1D* htempWalk1 = (TH1D*) gDirectory->FindObjectAny(
343 Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh));
344 if (NULL == htempWalk0 && NULL == htempWalk1) { // regenerate Walk histos
345 int iSide = 0;
346 htempWalk0 = new TH1D(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh),
347 Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]",
348 iSmType, iSm, iRpc, iCh, iSide),
349 nbClWalkBinX, fdTOTMin, fdTOTMax);
350 iSide = 1;
351 htempWalk1 = new TH1D(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh),
352 Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]",
353 iSmType, iSm, iRpc, iCh, iSide),
354 nbClWalkBinX, fdTOTMin, fdTOTMax);
355 }
356 if (NULL != htempWalk0 && NULL != htempWalk1) { // reinitialize Walk array
357 LOG(debug) << "Initialize Walk correction for "
358 << Form(" SmT%01d_sm%03d_rpc%03d_Ch%03d", iSmType, iSm, iRpc, iCh);
359 if (htempWalk0->GetNbinsX() != nbClWalkBinX)
360 LOG(error) << "CbmTaskTofClusterizerParWrite::InitCalibParameter: "
361 "Inconsistent Walk histograms";
362 for (int32_t iBin = 0; iBin < nbClWalkBinX; iBin++) {
363 vCPWalk[iCh][0][iBin] = htempWalk0->GetBinContent(iBin + 1);
364 vCPWalk[iCh][1][iBin] = htempWalk1->GetBinContent(iBin + 1);
365 if (iSmType == 0 && iSm == 0 && iRpc == 2 && iCh == 15) // debugging
366 LOG(info) << Form("Read new SmT%01d_sm%03d_rpc%03d_Ch%03d bin %d cen %f walk %f %f", iSmType, iSm,
367 iRpc, iCh, iBin, htempWalk0->GetBinCenter(iBin + 1), vCPWalk[iCh][0][iBin],
368 vCPWalk[iCh][1][iBin]);
369 if (5 == iSmType || 8 == iSmType) { // Pad structure, enforce consistency
370 if (vCPWalk[iCh][1][iBin] != vCPWalk[iCh][0][iBin]) {
371 LOG(fatal) << "Inconsisten walk values for TSRC " << iSmType << iSm << iRpc << iCh << ", Bin "
372 << iBin << ": " << vCPWalk[iCh][0][iBin] << ", " << vCPWalk[iCh][1][iBin];
373 }
374 vCPWalk[iCh][1][iBin] = vCPWalk[iCh][0][iBin];
375 htempWalk1->SetBinContent(iBin + 1, vCPWalk[iCh][1][iBin]);
376 }
377 }
378 }
379 else {
380 LOG(info) << "No Walk histograms for TSRC " << iSmType << iSm << iRpc << iCh;
381 }
382 }
383 // look for TcorY corrections
384 LOG(info) << "Check for TCorY in " << gDirectory->GetName();
385 TH1* hTCorY =
386 (TH1*) gDirectory->FindObjectAny(Form("calXY_SmT%d_sm%03d_rpc%03d_TOff_z_all_TcorY", iSmType, iSm, iRpc));
387 if (NULL != hTCorY) {
388 vCPTOffYBinWidth = hTCorY->GetBinWidth(0);
389 vCPTOffYRange = hTCorY->GetXaxis()->GetXmax();
390 LOG(info) << "Initialize TCorY: TSR " << iSmType << iSm << iRpc << ", Bwid " << vCPTOffYBinWidth
391 << ", Range " << vCPTOffYRange;
392 vCPTOffY.resize(hTCorY->GetNbinsX());
393 for (int iB = 0; iB < hTCorY->GetNbinsX(); iB++) {
394 vCPTOffY[iB] = hTCorY->GetBinContent(iB + 1);
395 }
396 }
397 }
398 else {
399 LOG(warning) << " Calibration histos " << Form("cl_SmT%01d_sm%03d_rpc%03d_XXX", iSmType, iSm, iRpc)
400 << " not found. ";
401 }
402 }
403 }
404
406 // <= To prevent histos from being sucked in by the param file of the TRootManager!
407 gFile = oldFile;
408 gDirectory = oldDir;
409 LOG(info) << "CbmTaskTofClusterizerParWrite::InitCalibParameter: initialization done";
410 return true;
411}
412
413bool CbmTaskTofClusterizerParWrite::InitAlgosTof()
414{
415 // Needed as external TOT values might be different than ones used for input histo reading (TO DO: FIX)
416 if (fTotMax != 0.) fdTOTMax = fTotMax;
417 if (fTotMin != 0.) fdTOTMin = fTotMin;
418 LOG(info) << "ToT init to Min " << fdTOTMin << " Max " << fdTOTMax;
419
421 gGeoManager->CdTop();
422
423 int32_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
424
425 // Create map with unique detector IDs and fill (needed only for dead strip array)
426 std::map<uint32_t, uint32_t> detIdIndexMap;
427 for (int32_t ind = 0; ind < fDigiBdfPar->GetNbDet(); ind++) {
428 int32_t iUniqueId = fDigiBdfPar->GetDetUId(ind);
429 detIdIndexMap[iUniqueId] = ind;
430 }
431
432 /* Hitfinding parameters */
433
434 cbm::algo::tof::HitfindSetup setup;
435 setup.rpcs.resize(iNbSmTypes);
436
437 for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
438
439 int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
440 int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
441 setup.NbSm.push_back(iNbSm);
442 setup.NbRpc.push_back(iNbRpc);
443
444 setup.rpcs[iSmType].resize(iNbSm * iNbRpc);
445
446 for (int32_t iSm = 0; iSm < iNbSm; iSm++) {
447 for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
448
449 cbm::algo::tof::HitfindSetup::Rpc par;
450
451 /* Cell geometry for hitfinding (can in principle be set channel-wise but is not needed for now) */
452
453 const int32_t rpcAddress = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
454 CbmTofCell* channelInfo = fDigiPar->GetCell(rpcAddress);
455 if (channelInfo == nullptr) {
456 continue;
457 }
458
459 // prepare local->global trafo
460 gGeoManager->FindNode(channelInfo->GetX(), channelInfo->GetY(), channelInfo->GetZ());
461 const double* tra_ptr = gGeoManager->MakePhysicalNode()->GetMatrix()->GetTranslation();
462 const double* rot_ptr = gGeoManager->GetCurrentMatrix()->GetRotationMatrix();
463 memcpy(par.cell.translation.data(), tra_ptr, 3 * sizeof(double));
464 memcpy(par.cell.rotation.data(), rot_ptr, 9 * sizeof(double));
465
466 par.cell.sizeX = channelInfo->GetSizex();
467 par.cell.sizeY = channelInfo->GetSizey();
468
469 /* Other hitfinding parameters */
470
471 const int32_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
472 const int32_t iDetIndx = detIdIndexMap[iDetId];
473 par.deadStrips = fvDeadStrips[iDetIndx];
474 par.posYMaxScal = fPosYMaxScal;
475 par.maxTimeDist = fdMaxTimeDist;
476 par.maxSpaceDist = fdMaxSpaceDist;
477 par.sigVel = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
478 par.timeRes = 0.08;
479 par.trackingStationId = fDigiBdfPar->GetTrackingStation(iSmType, iSm, iRpc);
480 par.CPTOffYBinWidth = fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc];
481 par.CPTOffY = fvCPTOffY[iSmType][iSm * iNbRpc + iRpc];
482 par.CPTOffYRange = fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc];
483
484 int32_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
485 par.chanPar.resize(iNbChan);
486
487 for (int32_t iCh = 0; iCh < iNbChan; iCh++) {
488 cbm::algo::tof::HitfindSetup::Channel& chan = par.chanPar[iCh];
489
490 const int32_t address = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
491 chan.address = address;
492 }
493 setup.rpcs[iSmType][iSm * iNbRpc + iRpc] = par;
494 }
495 }
496 }
497
498 /* Calibration parameters */
499
500 cbm::algo::tof::CalibrateSetup calib;
501 calib.rpcs.resize(iNbSmTypes);
502
503 for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
504 int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
505 int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
506 calib.NbSm.push_back(iNbSm);
507 calib.NbRpc.push_back(iNbRpc);
508 calib.rpcs[iSmType].resize(iNbSm * iNbRpc);
509
510 for (int32_t iSm = 0; iSm < iNbSm; iSm++) {
511 for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
512
513 cbm::algo::tof::CalibrateSetup::Rpc par;
514 par.numClWalkBinX = nbClWalkBinX;
515 par.TOTMax = fdTOTMax;
516 par.TOTMin = fdTOTMin;
517 par.swapChannelSides = fbSwapChannelSides;
518 par.channelDeadtime = fdChannelDeadtime;
519
520 int32_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
521 par.chanPar.resize(iNbChan);
522
523 for (int32_t iCh = 0; iCh < iNbChan; iCh++) {
524 cbm::algo::tof::CalibrateSetup::Channel& chan = par.chanPar[iCh];
525
526 chan.vCPTOff = fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh];
527 chan.vCPTotGain = fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh];
528 chan.vCPTotOff = fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh];
529 chan.vCPWalk = fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh];
530 }
531 calib.rpcs[iSmType][iSm * iNbRpc + iRpc] = par;
532 }
533 }
534 }
535
536 /* Write Yaml files */
537
538 {
539 namespace fs = boost::filesystem;
540 fs::path dirPath = fs::absolute(fs::weakly_canonical(fsRecoParOutputDir));
541
542 cbm::algo::yaml::Dump dump;
543 std::ofstream fout((dirPath / "TofHitfinderPar.yaml").string());
544 fout << dump(setup);
545 fout.close();
546
547 std::ofstream fcalout((dirPath / "TofCalibratePar.yaml").string());
548 fcalout << dump(calib);
549 fcalout.close();
550 }
551
552 return true;
553}
554
555
556bool CbmTaskTofClusterizerParWrite::InitAlgosBmon()
557{
558 // Needed as external TOT values might be different than ones used for input histo reading (TO DO: FIX)
559 if (fTotMax != 0.) fdTOTMax = fTotMax;
560 if (fTotMin != 0.) fdTOTMin = fTotMin;
561 LOG(info) << "ToT init to Min " << fdTOTMin << " Max " << fdTOTMax;
562
564 gGeoManager->CdTop();
565
566 // int32_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
567
568 // Create map with unique detector IDs and fill (needed only for dead strip array)
569 std::map<uint32_t, uint32_t> detIdIndexMap;
570 for (int32_t ind = 0; ind < fDigiBdfPar->GetNbDet(); ind++) {
571 int32_t iUniqueId = fDigiBdfPar->GetDetUId(ind);
572 detIdIndexMap[iUniqueId] = ind;
573 }
574
575 // Provide a selection mask between different diamonds
576 auto EstimateSelectionMask = [&](const auto& diamonds) -> uint32_t {
577 uint32_t selectionMask = 0;
578 for (auto itL = diamonds.begin(); itL != diamonds.end(); ++itL) {
579 for (auto itR = std::next(itL); itR != diamonds.end(); ++itR) {
580 selectionMask |= (itL->refAddress ^ itR->refAddress);
581 }
582 }
583 return selectionMask;
584 };
585
586
587 /* Hitfinding parameters */
588 int32_t iNbSm = fDigiBdfPar->GetNbSm(kBmonAssignedSmType);
589 int32_t iNbRpc = fDigiBdfPar->GetNbRpc(kBmonAssignedSmType);
590
591 cbm::algo::bmon::HitfindSetup setupHit;
592 cbm::algo::bmon::CalibrateSetup setupCal;
593 setupCal.diamonds.resize(iNbSm * iNbRpc);
594 setupHit.diamonds.resize(iNbSm * iNbRpc);
595
596 for (int32_t iSm = 0; iSm < iNbSm; iSm++) {
597 for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
598
599 const int32_t rpcAddress = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, kBmonAssignedSmType);
600 CbmTofCell* channelInfo = fDigiPar->GetCell(rpcAddress);
601 if (channelInfo == nullptr) {
602 continue;
603 }
604
605 //* Hitfinder
606 const int32_t iDetIndx = detIdIndexMap[rpcAddress];
607 cbm::algo::bmon::HitfindSetup::Diamond parHit;
608 parHit.refAddress = rpcAddress;
609 parHit.deadStrips = fvDeadStrips[iDetIndx];
610 parHit.maxTimeDist = fdMaxTimeDist;
611 parHit.timeRes = 0.08;
612
613 setupHit.diamonds[iSm * iNbRpc + iRpc] = parHit;
614
615 //* Calibration
616 cbm::algo::bmon::CalibrateSetup::Diamond parCal;
617 parCal.refAddress = rpcAddress;
618 parCal.numClWalkBinX = nbClWalkBinX;
619 parCal.TOTMax = fdTOTMax;
620 parCal.TOTMin = fdTOTMin;
621 parCal.channelDeadtime = fdChannelDeadtime;
622 int32_t iNbChan = fDigiBdfPar->GetNbChan(kBmonAssignedSmType, iRpc);
623 parCal.chanPar.resize(iNbChan);
624
625 for (int32_t iCh = 0; iCh < iNbChan; iCh++) {
626 cbm::algo::bmon::CalibrateSetup::Channel& chan = parCal.chanPar[iCh];
627
628 chan.vCPTOff = fvCPTOff[kBmonAssignedSmType][iSm * iNbRpc + iRpc][iCh][kBmonAssignedSide];
629 chan.vCPTotGain = fvCPTotGain[kBmonAssignedSmType][iSm * iNbRpc + iRpc][iCh][kBmonAssignedSide];
630 chan.vCPTotOff = fvCPTotOff[kBmonAssignedSmType][iSm * iNbRpc + iRpc][iCh][kBmonAssignedSide];
631 chan.vCPWalk = fvCPWalk[kBmonAssignedSmType][iSm * iNbRpc + iRpc][iCh][kBmonAssignedSide];
632 }
633
634 setupCal.diamonds[iSm * iNbRpc + iRpc] = parCal;
635 } // iRpc
636 } // iSm
637 setupHit.selectionMask = EstimateSelectionMask(setupHit.diamonds);
638 setupCal.selectionMask = setupHit.selectionMask;
639
640
641 /* Write Yaml files */
642
643 {
644 namespace fs = boost::filesystem;
645 fs::path dirPath = fs::absolute(fs::weakly_canonical(fsRecoParOutputDir));
646
647 cbm::algo::yaml::Dump dump;
648 std::ofstream fout((dirPath / "BmonHitfinderPar.yaml").string());
649 fout << dump(setupHit);
650 fout.close();
651
652 std::ofstream fcalout((dirPath / "BmonCalibratePar.yaml").string());
653 fcalout << dump(setupCal);
654 fcalout.close();
655 }
656
657 return true;
658}
hSvel
std::vector< TH1 * > hSvel
Definition CbmTofCalibrator.cxx:52
CbmTofDetectorId_v12b.h
CbmTofDetectorId_v14a.h
CbmTofDetectorId_v21a.h
k14a
@ k14a
Definition CbmTofGeoHandler.h:26
CalibrateSetup.h
Configuration of the calibrator for the BMON digis.
HitfindSetup.h
Parameters of the BMON hitfinder.
CbmTaskTofClusterizerParWrite::fvCPTOff
std::vector< std::vector< std::vector< std::vector< double > > > > fvCPTOff
Definition CbmTaskTofClusterizerParWrite.h:132
CbmTaskTofClusterizerParWrite::fvCPWalk
std::vector< std::vector< std::vector< std::vector< std::vector< double > > > > > fvCPWalk
Definition CbmTaskTofClusterizerParWrite.h:136
CbmTaskTofClusterizerParWrite::kBmonAssignedSide
static constexpr int kBmonAssignedSide
An RPC side, assigned to BMON diamonds.
Definition CbmTaskTofClusterizerParWrite.h:34
CbmTaskTofClusterizerParWrite::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition CbmTaskTofClusterizerParWrite.cxx:89
CbmTaskTofClusterizerParWrite::fvCPTOffYRange
std::vector< std::vector< double > > fvCPTOffYRange
Definition CbmTaskTofClusterizerParWrite.h:139
CbmTaskTofClusterizerParWrite::fvCPTOffY
std::vector< std::vector< std::vector< double > > > fvCPTOffY
Definition CbmTaskTofClusterizerParWrite.h:137
CbmTaskTofClusterizerParWrite::InitCalibParameter
bool InitCalibParameter()
Initialize other parameters not included in parameter classes.
Definition CbmTaskTofClusterizerParWrite.cxx:175
CbmTaskTofClusterizerParWrite::InitAlgosBmon
bool InitAlgosBmon()
Creates hit-finding and calibration parameters for BMON diamonds.
Definition CbmTaskTofClusterizerParWrite.cxx:556
CbmTaskTofClusterizerParWrite::InitParameters
bool InitParameters()
Initialize other parameters not included in parameter classes.
Definition CbmTaskTofClusterizerParWrite.cxx:113
CbmTaskTofClusterizerParWrite::InitAlgosTof
bool InitAlgosTof()
Create one algo object for each TOF RPC.
Definition CbmTaskTofClusterizerParWrite.cxx:413
CbmTaskTofClusterizerParWrite::kBmonAssignedSmType
static constexpr int kBmonAssignedSmType
A SmType, assigned to BMON diamonds.
Definition CbmTaskTofClusterizerParWrite.h:33
CbmTaskTofClusterizerParWrite::fvCPTOffYBinWidth
std::vector< std::vector< double > > fvCPTOffYBinWidth
Definition CbmTaskTofClusterizerParWrite.h:138
CbmTaskTofClusterizerParWrite::fvCPTotGain
std::vector< std::vector< std::vector< std::vector< double > > > > fvCPTotGain
Definition CbmTaskTofClusterizerParWrite.h:133
CbmTaskTofClusterizerParWrite::fvCPTotOff
std::vector< std::vector< std::vector< std::vector< double > > > > fvCPTotOff
Definition CbmTaskTofClusterizerParWrite.h:134
CbmTaskTofClusterizerParWrite::SetParContainers
virtual void SetParContainers()
Inherited from FairTask.
Definition CbmTaskTofClusterizerParWrite.cxx:99
CbmTofAddress::GetUniqueAddress
static uint32_t GetUniqueAddress(uint32_t Sm, uint32_t Rpc, uint32_t Channel, uint32_t Side=0, uint32_t SmType=0, uint32_t RpcType=0)
Definition CbmTofAddress.h:204
CbmTofCell::GetSizey
Double_t GetSizey() const
Definition CbmTofCell.h:40
CbmTofCell::GetY
Double_t GetY() const
Definition CbmTofCell.h:36
CbmTofCell::GetSizex
Double_t GetSizex() const
Definition CbmTofCell.h:39
CbmTofCell::GetX
Double_t GetX() const
Definition CbmTofCell.h:35
CbmTofCell::GetZ
Double_t GetZ() const
Definition CbmTofCell.h:37
CbmTofCreateDigiPar::Init
virtual InitStatus Init()
Definition CbmTofCreateDigiPar.cxx:113
CbmTofDetectorId_v14a
Definition CbmTofDetectorId_v14a.h:41
CbmTofDetectorId_v21a
Definition CbmTofDetectorId_v21a.h:41
CbmTofDigiBdfPar
Parameters class for the CBM ToF digitizer using beam data distributions.
Definition CbmTofDigiBdfPar.h:37
CbmTofGeoHandler::Init
Int_t Init(Bool_t isSimulation=kFALSE)
Definition CbmTofGeoHandler.cxx:46
cbm::algo::bmon::CalibrateSetup
BMON calibration per channel.
Definition bmon/CalibrateSetup.h:24
cbm::algo::bmon::HitfindSetup
Parameters for the BMON hitfinder.
Definition bmon/HitfindSetup.h:22
cbm::algo::bmon::HitfindSetup::diamonds
std::vector< Diamond > diamonds
Definition bmon/HitfindSetup.h:37
cbm::algo::tof::CalibrateSetup::Channel::vCPWalk
std::vector< std::vector< double > > vCPWalk
Definition tof/CalibrateSetup.h:26
cbm::algo::tof::CalibrateSetup::rpcs
std::vector< std::vector< Rpc > > rpcs
Definition tof/CalibrateSetup.h:53
cbm::algo::tof::HitfindSetup::Cell::translation
std::array< double, 3 > translation
Definition tof/HitfindSetup.h:27
cbm::algo::tof::HitfindSetup
Hitfind setup / Hardware cabling for TOF Used to create the hardware mapping for the TOF hitfinder.
Definition tof/HitfindSetup.h:22
cbm::algo::tof::HitfindSetup::rpcs
std::vector< std::vector< Rpc > > rpcs
Definition tof/HitfindSetup.h:72
cbm::algo::tof::HitfindSetup::NbSm
std::vector< int32_t > NbSm
Definition tof/HitfindSetup.h:70
cbm::algo::tof::HitfindSetup::NbRpc
std::vector< int32_t > NbRpc
Definition tof/HitfindSetup.h:71