CbmRoot
Loading...
Searching...
No Matches
CbmTrdRawToDigiR.cxx
Go to the documentation of this file.
1/* Copyright (C) 2020 Institut fuer Kernphysik, Goethe-Universitaet Frankfurt, Frankfurt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Etienne Bechtel [committer] */
4
5// Includes from TRD
6#include "CbmTrdRawToDigiR.h"
7
8#include "CbmTrdAddress.h"
9#include "CbmTrdDigi.h"
10
11#include "TMath.h"
12//#include "CbmSpadicRawMessage22.h"
13#include <TFile.h>
14#include <TProfile.h>
15#include <TProfile2D.h>
16
17#include <cassert>
18#include <chrono>
19#include <iostream>
20#include <string>
21#include <vector>
22
23using namespace std::chrono;
24
25//_________________________________________________________________________________
27 : TObject()
28 , fSampleMask()
30 , fElookupAsym()
31 , fElookupA()
32 , fElookupBig()
33 , fElookup()
34{
35}
36
37//_________________________________________________________________________________
39 : TObject()
40 , fSampleMask()
42 , fElookupAsym()
43 , fElookupA()
44 , fElookupBig()
45 , fElookup()
46{
47 SetReadFile(readfile);
48}
49
50
51//_________________________________________________________________________________
52CbmTrdRawToDigiR::CbmTrdRawToDigiR(Double_t cal, Double_t tau, Int_t mode)
53 : TObject()
54 , fSampleMask()
56 , fElookupAsym()
57 , fElookupA()
58 , fElookupBig()
59 , fElookup()
60{
61 SetCalibration(cal);
62 SetTau(tau);
63 SetRecoMode(mode);
64 SetPars(mode, cal, tau, fSampleMask);
65}
66
67
68//_________________________________________________________________________________
69CbmTrdRawToDigiR::CbmTrdRawToDigiR(Double_t cal, Double_t tau, std::vector<Int_t> mask)
70 : TObject()
71 , fSampleMask()
73 , fElookupAsym()
74 , fElookupA()
75 , fElookupBig()
76 , fElookup()
77{
78 SetCalibration(cal);
79 SetTau(tau);
80 SetRecoMask(mask);
81 SetPars(fRecoMode, cal, tau, mask);
82}
83
84
85void CbmTrdRawToDigiR::SetPars(Int_t mode, Double_t cal, Double_t tau, std::vector<Int_t> mask)
86{
87
88 //default
89 if (mode == 0) {
90 for (UInt_t i = 0; i <= mask.size(); i++)
91 fSampleMask.push_back(mask[i]);
92 fCalibration = cal;
93 fTau = tau;
94 Double_t sum = 0;
95 for (UInt_t i = 0; i < mask.size(); i++)
97 fEReco = sum;
99 }
100 if (mode == 1) {
101 for (Int_t i = fMaxBin; i <= fMaxBin; i++)
102 fSampleMask.push_back(i);
103 fCalibration = cal;
104 fTau = tau;
105 fLookUp = 0;
106 Double_t sum = 0;
107 for (UInt_t i = 0; i < fSampleMask.size(); i++)
109 fEReco = sum;
110 auto start = high_resolution_clock::now();
111 FillLookUps();
112 auto stop = high_resolution_clock::now();
113 auto duration = duration_cast<microseconds>(stop - start);
114 if (fDebug)
115 std::cout << std::endl
116 << std::endl
117 << "filling took: " << duration.count() << " microseconds " << std::endl
118 << fEReco << std::endl
119 << std::endl
120 << std::endl;
121 }
122 if (mode == 2) {
123 Double_t sum = 0;
124 for (UInt_t i = 0; i < fSampleMask.size(); i++)
126 fEReco = sum;
127 auto start = high_resolution_clock::now();
128 FillLookUps();
129 auto stop = high_resolution_clock::now();
130 auto duration = duration_cast<microseconds>(stop - start);
131 if (fDebug)
132 std::cout << std::endl
133 << std::endl
134 << "filling took: " << duration.count() << " microseconds " << std::endl
135 << fEReco << std::endl
136 << std::endl
137 << std::endl;
138 }
139}
140
141
143{
144 //default
145 Double_t sum = 0;
146 for (UInt_t i = 0; i < fSampleMask.size(); i++)
148 fEReco = sum;
149 auto start = high_resolution_clock::now();
150 if (fReadFile == "") FillLookUps(fWriteFile);
151 else
153 if (fLookUp == 4) FillLookUps("");
154 auto stop = high_resolution_clock::now();
155 auto duration = duration_cast<microseconds>(stop - start);
156 if (fDebug)
157 std::cout << std::endl
158 << std::endl
159 << "filling took: " << duration.count() << " microseconds " << std::endl
160 << fEReco << std::endl
161 << std::endl
162 << std::endl;
163 // if(fShapingOrder == 1) {SetMaxBin(2+fPresamples);SetMinBin(1+fPresamples);}
164 // if(fShapingOrder == 2) SetMaxBin(4+fPresamples);
165}
166
167
168void CbmTrdRawToDigiR::FillLookUps(std::string write)
169{
170
171 if (fLookUp == 1) {
172 for (Int_t shift = 0.; shift < CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC); shift++) {
173 Double_t sum = 0;
174 for (Int_t i = fMinBin; i <= fMaxBin; i++)
175 sum +=
177 fEReco = sum;
178 Float_t temp =
180
181 for (Int_t max = 0; max < fDynamicRange; max++) {
182 Float_t energy = max * 1.0 / temp;
183 Int_t a =
185 Int_t b = energy * fCalibration
187
188 Int_t mina = a - fExtrapolate * a;
189 Int_t maxa = a + fExtrapolate * a;
190 Int_t minb = b - fExtrapolate * b;
191 Int_t maxb = b + fExtrapolate * b;
192 if (!(fElookupAsym[max][a][b] > 0)) fElookupAsym[max][a][b] = shift;
193 fElookupSmall[shift][max] = energy;
194 if (fDebug) fQA->CreateHist("Asym", 512, -12.0, 500.0, 512, -12.0, 500.0);
195 if (fDebug && max == 200 && fQA->GetCont2D("Asym", a, b) == 0.) fQA->Fill("Asym", a, b, shift);
196 if (fDebug) fQA->CreateHist("Look", 63, -0.5, 62.5, 512, -12.0, 500.0);
197 if (fDebug) fQA->Fill("Look", shift, max, energy);
198 Int_t extrapolate = 0;
199 while (true) {
200 extrapolate++;
201 if ((b - extrapolate) <= minb || (a - extrapolate) <= mina) break;
202 if (fElookupAsym[max][a - extrapolate][b - extrapolate] > 0) continue;
203 fElookupAsym[max][a - extrapolate][b - extrapolate] = fElookupAsym[max][a][b];
204 if (fDebug && max == 200) fQA->Fill("Asym", a - extrapolate, b - extrapolate, fElookupAsym[max][a][b]);
205 Int_t count = 1;
206 if (!(fElookupAsym[max][a - extrapolate - count][b - extrapolate] > 0)) {
207 fElookupAsym[max][a - extrapolate - count][b - extrapolate] =
208 fElookupAsym[max][a - extrapolate][b - extrapolate];
209 if (fDebug && max == 200)
210 fQA->Fill("Asym", a - extrapolate - count, b - extrapolate, fElookupAsym[max][a][b]);
211 count++;
212 }
213 count = 1;
214 if (!(fElookupAsym[max][a - extrapolate + count][b - extrapolate] > 0)) {
215 fElookupAsym[max][a - extrapolate + count][b - extrapolate] =
216 fElookupAsym[max][a - extrapolate][b - extrapolate];
217 if (fDebug && max == 200)
218 fQA->Fill("Asym", a - extrapolate + count, b - extrapolate, fElookupAsym[max][a][b]);
219 count++;
220 }
221 }
222 extrapolate = 0;
223 while (true) {
224 extrapolate++;
225 if ((b + extrapolate) >= maxb || (a + extrapolate) >= maxa) break;
226 if (fElookupAsym[max][a + extrapolate][b + extrapolate] > 0) continue;
227 fElookupAsym[max][a + extrapolate][b + extrapolate] = fElookupAsym[max][a][b];
228 if (fDebug && max == 200) fQA->Fill("Asym", a + extrapolate, b + extrapolate, fElookupAsym[max][a][b]);
229 Int_t count = 1;
230 if (!(fElookupAsym[max][a + extrapolate - count][b + extrapolate] > 0)) {
231 fElookupAsym[max][a + extrapolate - count][b + extrapolate] =
232 fElookupAsym[max][a + extrapolate][b + extrapolate];
233 if (fDebug && max == 200)
234 fQA->Fill("Asym", a + extrapolate - count, b + extrapolate, fElookupAsym[max][a][b]);
235 count++;
236 }
237 count = 1;
238 if (!(fElookupAsym[max][a + extrapolate + count][b + extrapolate] > 0)) {
239 fElookupAsym[max][a + extrapolate + count][b + extrapolate] =
240 fElookupAsym[max][a + extrapolate][b + extrapolate];
241 if (fDebug && max == 200)
242 fQA->Fill("Asym", a + extrapolate + count, b + extrapolate, fElookupAsym[max][a][b]);
243 count++;
244 }
245 }
246 }
247 }
248 if (write != "") WriteMaps(write);
249 }
250 if (fLookUp == 2) {
251 for (Int_t shift = 0.; shift < CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC); shift++) {
252 Double_t sum = 0;
253 for (Int_t i = fMinBin; i <= fMaxBin; i++)
254 sum +=
256 fEReco = sum;
257 Float_t temp =
259 Int_t extrapolate = 0;
260 for (Int_t max = 0; max < fDynamicRange; max++) {
261 Float_t energy = max * 1.0 / temp;
262 Int_t a =
264 Int_t b = energy * fCalibration
266
267 Int_t mina = a - fExtrapolate * a;
268 Int_t maxa = a + fExtrapolate * a;
269 Int_t minb = b - fExtrapolate * b;
270 Int_t maxb = b + fExtrapolate * b;
271 if (!(fElookupAsym[max][a][b] > 0)) fElookupAsym[max][a][b] = shift;
272 sum = 0.;
273 for (UInt_t i = 0; i < fSampleMask.size(); i++)
274 sum += energy * fCalibration
276
277 fElookupSmall[shift][sum] = energy;
278 if (fDebug) fQA->CreateHist("Asym", 512, -12.0, 500.0, 512, -12.0, 500.0);
279 if (fDebug && max == 200 && fQA->GetCont2D("Asym", a, b) == 0.) fQA->Fill("Asym", a, b, shift);
280 if (fDebug) fQA->CreateHist("Look", 63, -0.5, 62.5, 512, -12.0, 2000.0);
281 if (fDebug) fQA->Fill("Look", shift, sum, energy);
282 extrapolate = 0;
283 while (true) {
284 extrapolate++;
285 if ((b - extrapolate) <= minb || (a - extrapolate) <= mina) break;
286 if (fElookupAsym[max][a - extrapolate][b - extrapolate] > 0) continue;
287 fElookupAsym[max][a - extrapolate][b - extrapolate] = fElookupAsym[max][a][b];
288 if (fDebug && max == 200) fQA->Fill("Asym", a - extrapolate, b - extrapolate, fElookupAsym[max][a][b]);
289 Int_t count = 1;
290 if (!(fElookupAsym[max][a - extrapolate - count][b - extrapolate] > 0)) {
291 fElookupAsym[max][a - extrapolate - count][b - extrapolate] =
292 fElookupAsym[max][a - extrapolate][b - extrapolate];
293 if (fDebug && max == 200)
294 fQA->Fill("Asym", a - extrapolate - count, b - extrapolate, fElookupAsym[max][a][b]);
295 count++;
296 }
297 count = 1;
298 if (!(fElookupAsym[max][a - extrapolate + count][b - extrapolate] > 0)) {
299 fElookupAsym[max][a - extrapolate + count][b - extrapolate] =
300 fElookupAsym[max][a - extrapolate][b - extrapolate];
301 if (fDebug && max == 200)
302 fQA->Fill("Asym", a - extrapolate + count, b - extrapolate, fElookupAsym[max][a][b]);
303 count++;
304 }
305 }
306 extrapolate = 0;
307 while (true) {
308 extrapolate++;
309 if ((b + extrapolate) >= maxb || (a + extrapolate) >= maxa) break;
310 if (fElookupAsym[max][a + extrapolate][b + extrapolate] > 0) continue;
311 fElookupAsym[max][a + extrapolate][b + extrapolate] = fElookupAsym[max][a][b];
312 if (fDebug && max == 200) fQA->Fill("Asym", a + extrapolate, b + extrapolate, fElookupAsym[max][a][b]);
313 Int_t count = 1;
314 if (!(fElookupAsym[max][a + extrapolate - count][b + extrapolate] > 0)) {
315 fElookupAsym[max][a + extrapolate - count][b + extrapolate] =
316 fElookupAsym[max][a + extrapolate][b + extrapolate];
317 if (fDebug && max == 200)
318 fQA->Fill("Asym", a + extrapolate - count, b + extrapolate, fElookupAsym[max][a][b]);
319 count++;
320 }
321 count = 1;
322 if (!(fElookupAsym[max][a + extrapolate + count][b + extrapolate] > 0)) {
323 fElookupAsym[max][a + extrapolate + count][b + extrapolate] =
324 fElookupAsym[max][a + extrapolate][b + extrapolate];
325 if (fDebug && max == 200)
326 fQA->Fill("Asym", a + extrapolate + count, b + extrapolate, fElookupAsym[max][a][b]);
327 count++;
328 }
329 }
330 }
331 }
332 if (write != "") WriteMaps(write);
333 }
334 if (fLookUp == 4) {
335 for (Int_t max = 0; max < fDynamicRange; max++) {
336 Float_t energy =
338 fElookup[max] = energy;
339 }
340 }
341 else {
342 Int_t range = fDynamicRange * fSampleMask.size();
343 for (Int_t n = 0; n <= range; n++) {
344 Float_t value = n / fEReco;
345 fElookup[n] = value;
346 }
347 }
348}
349
350
351void CbmTrdRawToDigiR::ReadMaps(std::string file)
352{
353 if (fLookUp == 3) {
354 TFile f(file.data(), "OPEN");
355 LOG_IF(fatal, !f.IsOpen()) << "parameter file " << file << " does not exist!";
356 TProfile2D* h;
357 h = f.Get<TProfile2D>("MAX ADC");
358 LOG_IF(fatal, !h) << "No histogram MAX ADC found in file " << file;
359 for (Int_t x = 1; x <= h->GetNbinsX(); x++) {
360 for (Int_t y = 1; y <= h->GetNbinsY(); y++) {
361 fElookupSmall[h->GetXaxis()->GetBinCenter(x)][h->GetYaxis()->GetBinCenter(y)] = h->GetBinContent(x, y);
362 if (fDebug) fQA->CreateHist("Read Small", 63, -0.5, 62.5, 512, -12., 500.);
363 if (fDebug)
364 fQA->Fill("Read Small", h->GetXaxis()->GetBinCenter(x), h->GetYaxis()->GetBinCenter(y),
365 h->GetBinContent(x, y));
366 }
367 }
368 h = f.Get<TProfile2D>("ASYM MAP");
369 LOG_IF(fatal, !h) << "No histogram ASYM MAP found in file " << file;
370 for (Int_t x = 1; x <= h->GetNbinsX(); x++) {
371 for (Int_t y = 1; y <= h->GetNbinsY(); y++) {
372 fElookupA[h->GetXaxis()->GetBinCenter(x)][h->GetYaxis()->GetBinCenter(y)] = h->GetBinContent(x, y);
373 if (fDebug) fQA->CreateHist("Read Asym", 512, 0., 512., 512, -12., 500.);
374 if (fDebug)
375 fQA->Fill("Read Asym", h->GetXaxis()->GetBinCenter(x), h->GetYaxis()->GetBinCenter(y),
376 h->GetBinContent(x, y));
377 }
378 }
379 f.Close();
380 }
381 if (fLookUp == 4) {
382 TFile f(file.data(), "OPEN");
383 LOG_IF(fatal, !f.IsOpen()) << "parameter file " << file << " does not exist!";
384 TProfile2D* h;
385 h = f.Get<TProfile2D>("MAX ADC");
386 LOG_IF(fatal, !h) << "No histogram MAX ADC found in file " << file;
387 for (Int_t x = 1; x <= h->GetNbinsX(); x++) {
388 for (Int_t y = 1; y <= h->GetNbinsY(); y++) {
389 fElookupSmall[h->GetXaxis()->GetBinCenter(x)][h->GetYaxis()->GetBinCenter(y)] = h->GetBinContent(x, y);
390 if (fDebug) fQA->CreateHist("Read Small", 63, -0.5, 62.5, 512, -12., 500.);
391 if (fDebug)
392 fQA->Fill("Read Small", h->GetXaxis()->GetBinCenter(x), h->GetYaxis()->GetBinCenter(y),
393 h->GetBinContent(x, y));
394 }
395 }
396 h = f.Get<TProfile2D>("ASYM MAP");
397 LOG_IF(fatal, !h) << "No histogram ASYM MAP found in file " << file;
398 for (Int_t x = 1; x <= h->GetNbinsX(); x++) {
399 for (Int_t y = 1; y <= h->GetNbinsY(); y++) {
400 fElookupA[h->GetXaxis()->GetBinCenter(x)][h->GetYaxis()->GetBinCenter(y)] = h->GetBinContent(x, y);
401 if (fDebug) fQA->CreateHist("Read Asym", 512, 0., 512., 512, -12., 500.);
402 if (fDebug)
403 fQA->Fill("Read Asym", h->GetXaxis()->GetBinCenter(x), h->GetYaxis()->GetBinCenter(y),
404 h->GetBinContent(x, y));
405 }
406 }
407 f.Close();
408 }
409}
410
411CbmTrdDigi* CbmTrdRawToDigiR::MakeDigi(std::vector<Int_t> samples, Int_t channel, Int_t uniqueModuleId, ULong64_t time,
412 Bool_t FN)
413{
414 Float_t digicharge = 0;
415 Int_t samplesum = 0;
416 for (size_t i = 0; i < fSampleMask.size(); i++) {
417 samplesum += samples[fSampleMask[i]];
418 }
419 if (fReadFile == "") {
420 if (fLookUp == 1) {
421 digicharge = fElookupSmall[fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]][samples[fMaxBin]];
422 }
423 if (fLookUp == 2) {
424 digicharge = fElookupSmall[fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]][samplesum];
425 }
426 else {
427 digicharge = fElookup[samplesum];
428 }
429 }
430 else {
431 if (fLookUp == 3) { digicharge = fElookupSmall[fElookupA[samples[fMinBin]][samples[fMaxBin]]][samples[fMaxBin]]; }
432 if (fLookUp == 4 && !FN) {
433 digicharge = fElookupSmall[fElookupA[samples[fMinBin]][samples[fMaxBin]]][samples[fMaxBin]];
434 }
435 if (fLookUp == 4 && FN) { digicharge = fElookup[samples[fMaxBin]]; }
436 }
437 // The triggertype is set later by the class (moduleSimR) using this function
438 CbmTrdDigi* digi = new CbmTrdDigi(channel, uniqueModuleId, digicharge, time, CbmTrdDigi::eTriggerType::kNTrg, 0);
439
440 return digi;
441}
442
443Float_t CbmTrdRawToDigiR::GetTimeShift(std::vector<Int_t> samples)
444{
445
446 Float_t shift = 0;
447 if (fReadFile == "") {
448 if (fLookUp == 1) { shift = fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]; }
449 else {
450 return 0.;
451 }
452 }
453 else {
454 if (fLookUp == 3) { shift = fElookupA[samples[fMinBin]][samples[fMaxBin]]; }
455 if (fLookUp == 4) { shift = fElookupA[samples[fMinBin]][samples[fMaxBin]]; }
456 }
457
458 return shift;
459}
460
461Double_t CbmTrdRawToDigiR::GetCharge(std::vector<Int_t> samples, Int_t shift)
462{
463
464 Float_t charge = 0;
465 if (fReadFile == "") {
466 if (fLookUp == 1) {
467 if (shift > -1) charge = fElookupSmall[shift][samples[fMaxBin]];
468 else
469 charge = fElookupSmall[fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]][samples[fMaxBin]];
470 }
471 if (fLookUp == 2) {
472 Int_t samplesum = 0;
473 for (size_t i = 0; i < fSampleMask.size(); i++) {
474 samplesum += samples[fSampleMask[i]];
475 }
476 if (shift > -1) charge = fElookupSmall[shift][samplesum];
477 else
478 charge = fElookupSmall[fElookupAsym[samples[fMaxBin]][samples[fMinBin]][samples[fHighBin]]][samplesum];
479 }
480 else {
481 return 0.;
482 }
483 }
484 else {
485 if (fLookUp == 3) {
486 if (shift > -1) charge = fElookupSmall[shift][samples[fMaxBin]];
487 else
488 charge = fElookupSmall[fElookupA[samples[fMinBin]][samples[fMaxBin]]][samples[fMaxBin]];
489 }
490 if (fLookUp == 4) {
491 if (shift > -1) charge = fElookupSmall[shift][samples[fMaxBin]];
492 else
493 charge = fElookupSmall[fElookupA[samples[fMinBin]][samples[fMaxBin]]][samples[fMaxBin]];
494 }
495 }
496
497 return charge;
498}
499
500
501//_______________________________________________________________________________
503{
504
505 if (fShapingOrder == 1) return (t / fTau) * TMath::Exp(-(t / fTau));
506 if (fShapingOrder == 2) return (t / fTau) * (t / fTau) * TMath::Exp(-(t / fTau));
507
508 return 0.;
509}
Helper class to convert unique channel ID back and forth.
friend fscal max(fscal x, fscal y)
float Float_t
int Int_t
bool Bool_t
Generates beam ions for transport simulation.
static float Clk(eCbmTrdAsicType ty)
DAQ clock accessor for each ASIC.
Definition CbmTrdDigi.h:110
void FillLookUps(std::string write="")
CbmTrdDigi * MakeDigi(std::vector< Int_t > samples, Int_t channel, Int_t uniqueModuleId, ULong64_t time, Bool_t FN=false)
Double_t CalcResponse(Double_t t)
void SetPars(Int_t mode, Double_t cal, Double_t tau, std::vector< Int_t > mask)
std::map< Int_t, std::map< Int_t, std::map< Int_t, Int_t > > > fElookupBig
std::map< Int_t, std::map< Int_t, Float_t > > fElookupSmall
std::map< Int_t, Float_t > fElookup
void SetReadFile(std::string file)
std::map< Int_t, std::map< Int_t, Int_t > > fElookupA
void SetRecoMask(std::vector< Int_t > mask)
void WriteMaps(std::string file="")
CbmTrdRawToDigiR()
default Constructor with messages
std::map< Int_t, std::map< Int_t, std::map< Int_t, Int_t > > > fElookupAsym
void SetCalibration(Double_t cal)
std::string fReadFile
void ReadMaps(std::string file="")
Double_t GetCharge(std::vector< Int_t > samples, Int_t shift=-1)
void SetTau(Double_t tau)
std::string fWriteFile
void SetRecoMode(Int_t mode)
Float_t GetTimeShift(std::vector< Int_t > samples)
std::vector< Int_t > fSampleMask
CbmTrdCheckUtil * fQA