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