CbmRoot
Loading...
Searching...
No Matches
CbmTrdModuleSimR.cxx
Go to the documentation of this file.
1/* Copyright (C) 2018-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Etienne Bechtel, Florian Uhlig [committer] */
4
5// Includes from TRD
6#include "CbmTrdModuleSimR.h"
7
8#include "CbmTrdAddress.h"
9#include "CbmTrdDigi.h"
10#include "CbmTrdDigitizer.h"
11#include "CbmTrdParModAsic.h"
12#include "CbmTrdParModDigi.h"
13#include "CbmTrdParSpadic.h"
14#include "CbmTrdPoint.h"
15#include "CbmTrdRadiator.h"
16
17// Includes from CbmRoot
18//#include "CbmDaqBuffer.h"
19#include "CbmMatch.h"
20
21// Includes from FairRoot
22#include <Logger.h>
23
24// Includes from Root
25#include <TClonesArray.h>
26#include <TGeoManager.h>
27#include <TMath.h>
28#include <TRandom3.h>
29#include <TStopwatch.h>
30#include <TVector3.h>
31
32// Includes from C++
33#include "CbmTrdRawToDigiR.h"
34
35#include <cmath>
36#include <fstream>
37#include <iomanip>
38#include <iostream>
39
40
41//_________________________________________________________________________________
43 : CbmTrdModuleSim(mod, ly, rot)
44 , fSigma_noise_keV(0.2)
45 , fMinimumChargeTH(.5e-06)
46 , fCurrentTime(-1.)
47 , fAddress(-1.)
48 , fLastEventTime(-1)
49 , fEventTime(-1)
50 , fLastTime(-1)
51 , fCollectTime(2048)
52 , //in pulsed mode
55 ,
56
57 fPulseSwitch(true)
58 , fPrintPulse(false)
59 , //only for debugging
60 fTimeShift(true)
61 , //for pulsed mode
62 fAddCrosstalk(true)
63 , //for pulsed mode
64 fClipping(true)
65 , //for pulsed mode
66 fepoints(0)
67 , fAdcNoise(2)
69 , fCrosstalkLevel(0.01)
70 ,
71
75 ,
76
78 , fPulseBuffer()
79 , fMultiBuffer()
80 , fTimeBuffer()
81 , fShiftQA()
82 , fLinkQA()
83 , fMCQA()
84 , fMCBuffer()
85
86{
87 FairRootManager* ioman = FairRootManager::Instance();
88 fTimeSlice = (CbmTimeSlice*) ioman->GetObject("TimeSlice.");
89
91 SetNameTitle(Form("TrdSimR%d", mod), "Simulator for rectangular read-out.");
92
94 TFile* oldFile = gFile;
95 TDirectory* oldDir = gDirectory;
96
97 TString dir = getenv("VMCWORKDIR");
98 TString filename = dir + "/parameters/trd/FeatureExtractionLookup.root";
99 TFile* f = new TFile(filename, "OPEN");
100 LOG_IF(fatal, !f->IsOpen()) << "parameter file " << filename << " does not exist!";
101 fDriftTime = f->Get<TH2D>("Drift");
102 LOG_IF(fatal, !fDriftTime) << "No histogram Drift founfd in file " << filename;
103 fDriftTime->SetDirectory(0);
104 f->Close();
105
107 gFile = oldFile;
108 gDirectory = oldDir;
109
110 oldFile = nullptr;
111 delete oldFile;
112
113 if (fPulseSwitch) {
116 }
117}
118
119//_______________________________________________________________________________
120void CbmTrdModuleSimR::AddDigi(Int_t address, Double_t charge, Double_t /*chargeTR*/, Double_t time, Int_t trigger)
121{
122 Double_t weighting = charge;
124 TVector3 padPos, padPosErr;
125 fDigiPar->GetPadPosition(address, padPos, padPosErr);
126 Double_t distance = sqrt(pow(fXYZ[0] - padPos[0], 2) + pow(fXYZ[1] - padPos[1], 2));
127 weighting = 1. / distance;
128 }
129
130 // std::map<Int_t, std::pair<CbmTrdDigi*, CbmMatch*> >::iterator it = fDigiMap.find(address);
131 // std::map<Int_t, std::pair<CbmTrdDigi*, CbmMatch*> >::iterator it;
132
133 Int_t col = CbmTrdAddress::GetColumnId(address);
134 Int_t row = CbmTrdAddress::GetRowId(address);
135 Int_t sec = CbmTrdAddress::GetSectorId(address);
136 Int_t ncols = fDigiPar->GetNofColumns();
137 Int_t channel(0);
138 for (Int_t isec(0); isec < sec; isec++)
139 channel += ncols * fDigiPar->GetNofRowsInSector(isec);
140 channel += ncols * row + col;
141
142 std::map<Int_t, std::pair<CbmTrdDigi*, CbmMatch*>>::iterator it = fDigiMap.find(address);
143 if (it == fDigiMap.end()) { // Pixel not yet in map -> Add new pixel
144 CbmMatch* digiMatch = new CbmMatch();
145 digiMatch->AddLink(CbmLink(weighting, fPointId, fEventId, fInputId));
146 AddNoise(charge);
147
149 if (trigger == 1) triggertype = CbmTrdDigi::eTriggerType::kSelf;
150 if (trigger == 2) triggertype = CbmTrdDigi::eTriggerType::kNeighbor;
151
152 CbmTrdDigi* digi = new CbmTrdDigi(channel, fModAddress, charge * 1e6, ULong64_t(time), triggertype, 0);
153
154 digi->SetFlag(0, true);
155 if (fDigiPar->GetPadSizeY(sec) == 1.5) digi->SetErrorClass(1);
156 if (fDigiPar->GetPadSizeY(sec) == 4.) digi->SetErrorClass(2);
157 if (fDigiPar->GetPadSizeY(sec) == 6.75) digi->SetErrorClass(3);
158 if (fDigiPar->GetPadSizeY(sec) == 12.) digi->SetErrorClass(4);
159
160 fDigiMap[address] = std::make_pair(digi, digiMatch);
161
162 it = fDigiMap.find(address);
163 // it->second.first->SetAddressModule(fModAddress);//module); <- now handled in the digi contructor
164 }
165 else {
166 it->second.first->AddCharge(charge * 1e6);
167 it->second.second->AddLink(CbmLink(weighting, fPointId, fEventId, fInputId));
168 }
169}
170
171//_______________________________________________________________________________
173{
174
175 std::vector<std::pair<CbmTrdDigi*, CbmMatch*>> analog = fAnalogBuffer[address];
176 std::vector<std::pair<CbmTrdDigi*, CbmMatch*>>::iterator it;
177 CbmTrdDigi* digi = analog.begin()->first;
178 CbmMatch* digiMatch = new CbmMatch();
179 //printf("CbmTrdModuleSimR::ProcessBuffer(%10d)=%3d\n", address, digi->GetAddressChannel());
180
181 // Int_t trigger = 0;
182 Float_t digicharge = 0;
183 // Float_t digiTRcharge=0;
184 for (it = analog.begin(); it != analog.end(); it++) {
185 digicharge += it->first->GetCharge();
186 digiMatch->AddLink(it->second->GetLink(0));
187 //printf(" add charge[%f] trigger[%d]\n", it->first->GetCharge(), it->first->GetTriggerType());
188 }
189 digi->SetCharge(digicharge);
190 digi->SetTriggerType(fAnalogBuffer[address][0].first->GetTriggerType());
191
192 // std::cout<<digicharge<<std::endl;
193
194 // if(analog.size()>1) for (it=analog.begin()+1 ; it != analog.end(); it++) if(it->first) delete it->first;
195
196 fDigiMap[address] = std::make_pair(digi, digiMatch);
197
198 fAnalogBuffer.erase(address);
199}
200
201//_______________________________________________________________________________
202void CbmTrdModuleSimR::ProcessPulseBuffer(Int_t address, Bool_t FNcall, Bool_t MultiCall = false, Bool_t down = true,
203 Bool_t up = true)
204{
205
206 std::map<Int_t, std::pair<std::vector<Double_t>, CbmMatch*>>::iterator iBuff = fPulseBuffer.find(address);
207 std::map<Int_t, Double_t>::iterator tBuff = fTimeBuffer.find(address);
208
209 if (iBuff == fPulseBuffer.end() || tBuff == fTimeBuffer.end()) return;
210 Int_t trigger = CheckTrigger(fPulseBuffer[address].first);
211 if (fPulseBuffer[address].first.size() < 32) {
212 return;
213 }
214
215 if (trigger == 0 && !FNcall) {
216 return;
217 }
218 if (trigger == 1 && FNcall) {
219 FNcall = false;
220 }
221
222 Int_t col = CbmTrdAddress::GetColumnId(address);
223 Int_t row = CbmTrdAddress::GetRowId(address);
224 Int_t sec = CbmTrdAddress::GetSectorId(address);
225 Int_t ncols = fDigiPar->GetNofColumns();
226 Int_t channel(0);
227 for (Int_t isec(0); isec < sec; isec++)
228 channel += ncols * fDigiPar->GetNofRowsInSector(isec);
229 channel += ncols * row + col;
230
231
232 CbmMatch* digiMatch = new CbmMatch();
233 std::vector<Int_t> temp;
234 for (size_t i = 0; i < fPulseBuffer[address].first.size(); i++) {
235 Int_t noise = AddNoiseADC();
236 Int_t cross = AddCrosstalk(address, i, sec, row, col, ncols);
237
238 fPulseBuffer[address].first[i] = fPulseBuffer[address].first[i] + noise + cross;
239 if (fPulseBuffer[address].first[i] <= fClipLevel && fPulseBuffer[address].first[i] >= -12.)
240 temp.push_back((Int_t) fPulseBuffer[address].first[i]);
241 if (fPulseBuffer[address].first[i] > fClipLevel) temp.push_back(fClipLevel - 1);
242 if (fPulseBuffer[address].first[i] < -12.) temp.push_back(-12.);
243 if (fDebug) fQA->Fill("Pulse", i, fPulseBuffer[address].first[i]);
244 if (fDebug) fQA->CreateHist("Pulse self", 32, -0.5, 31.5, 512, -12., 500.);
245 if (fDebug) fQA->CreateHist("Pulse FN", 32, -0.5, 31.5, 512, -12., 500.);
246 if (fDebug && trigger == 1) fQA->Fill("Pulse self", i, fPulseBuffer[address].first[i]);
247 if (fDebug && trigger == 0) fQA->Fill("Pulse FN", i, fPulseBuffer[address].first[i]);
248 }
249
250 Float_t shift = fMessageConverter->GetTimeShift(temp);
251 Float_t corr = fMinDrift; //correction of average sampling to clock difference and systematic average drift time
253 //correction for EB case is done later, due to the event time 0 and the unsigned data member for the time in the digi
254
255 if (fTimeSlice) {
256 if (fTimeBuffer[address] + corr - shift < fTimeSlice->GetStartTime()) {
257 delete digiMatch;
258 delete fPulseBuffer[address].second;
259 fPulseBuffer.erase(address);
260 fTimeBuffer.erase(address);
261 return;
262 }
263 }
264
265 CbmTrdDigi* digi = nullptr;
266 if (fTimeBuffer[address] + corr - shift > 0.)
267 digi = fMessageConverter->MakeDigi(temp, channel, fModAddress, fTimeBuffer[address] + corr - shift, true);
268 else
269 digi = fMessageConverter->MakeDigi(temp, channel, fModAddress, fTimeBuffer[address] + corr, true);
270
271 if (fDigiPar->GetPadSizeY(sec) == 1.5) digi->SetErrorClass(1);
272 if (fDigiPar->GetPadSizeY(sec) == 4.) digi->SetErrorClass(2);
273 if (fDigiPar->GetPadSizeY(sec) == 6.75) digi->SetErrorClass(3);
274 if (fDigiPar->GetPadSizeY(sec) == 12.) digi->SetErrorClass(4);
275 if (!CbmTrdDigitizer::IsTimeBased()) digi->SetFlag(1, true);
276
277 Int_t shiftcut = fShiftQA[address] * 10;
278 Float_t timeshift = shiftcut / 10.0;
279 if (temp[fMinBin] == fClipLevel - 1 && temp[fMaxBin] == fClipLevel - 1) digi->SetCharge(35.);
280
281 std::vector<CbmLink> links = fPulseBuffer[address].second->GetLinks();
282 Int_t Links = 0.;
283 for (UInt_t iLinks = 0; iLinks < links.size(); iLinks++) {
284 digiMatch->AddLink(links[iLinks]);
285 Links++;
286 }
287
288 if (fDebug) {
289 fQA->CreateHist("MC", 200, 0., 50.);
290 fQA->Fill("MC", fMCQA[address]);
291 fQA->CreateHist("rec shift", 63, -0.5, 62.5);
292 fQA->Fill("rec shift", shift);
293 fQA->CreateHist("T res", 70, -35, 35);
294 fQA->Fill("T res", shift - timeshift);
295 fQA->CreateHist("Time", 100, -50, 50);
296 fQA->Fill("Time", digi->GetTime() - fLinkQA[address][0]["Time"]);
297
298
299 fQA->CreateProfile("MAX ADC", 63, -0.5, 62.5, 512, -12., 500.);
300 fQA->FillProfile("MAX ADC", timeshift, temp[fMaxBin], fMCQA[address]);
301 fQA->CreateProfile("ASYM MAP", 512, -12., 500., 512, -12., 500.);
302 fQA->FillProfile("ASYM MAP", temp[fMinBin], temp[fMaxBin], timeshift);
303
304
305 if (trigger == 1 && MultiCall && digi->GetCharge() > 0.)
306 fQA->Fill("Multi Quote", 1);
307 else
308 fQA->Fill("Multi Quote", 0);
309
310 if (trigger == 1 && !MultiCall) {
311 fQA->CreateHist("E self", 200, 0., 50.0);
312 fQA->Fill("E self", digi->GetCharge());
313 fQA->CreateHist("E res", 400., -1., 1.);
314 fQA->Fill("E res", digi->GetCharge() - fMCQA[address]);
315 fQA->CreateHist("E res rel", 400., -1., 1.);
316 fQA->Fill("E res rel", (digi->GetCharge() - fMCQA[address]) / fMCQA[address]);
317 fQA->CreateHist("Charge Max", 200, 0., 50.0, 512, -12., 500.);
318 fQA->Fill("Charge Max", fMessageConverter->GetCharge(temp, timeshift), temp[fMaxBin]);
319 fQA->CreateHist("Links Res Self", 400., -1., 1., 5, -0.5, 4.5);
320 fQA->Fill("Links Res Self", (digi->GetCharge() - fMCQA[address]) / fMCQA[address], digiMatch->GetNofLinks());
321
322 if (Links == 1) {
323 fQA->CreateProfile("1 Link Prof", 32, 0., 32.);
324 for (size_t i = 0; i < fPulseBuffer[address].first.size(); i++) {
325 fQA->FillProfile("1 Link Prof", i, temp[i]);
326 }
327 }
328
329 if (fLinkQA[address].size() == 2 && fLinkQA[address][1]["Event"] != fLinkQA[address][0]["Event"]) {
330 fQA->CreateHist("Links QA time", 400., -1., 1., 500, 0., 2000.);
331 fQA->Fill("Links QA time", (digi->GetCharge() - fMCQA[address]) / fMCQA[address],
332 fLinkQA[address][1]["Time"] - fLinkQA[address][0]["Time"]);
333 fQA->CreateProfile("2 Link Prof", 32, 0., 32., 100, 0., 2000.);
334 for (size_t i = 0; i < fPulseBuffer[address].first.size(); i++) {
335 fQA->FillProfile("2 Link Prof", i, fLinkQA[address][1]["Time"] - fLinkQA[address][0]["Time"], temp[i]);
336 }
337 }
338
339 if (Links == 2 && links[0].GetEntry() == links[1].GetEntry()) {
340 fQA->CreateHist("In Event Res", 400., -1., 1.);
341 fQA->Fill("In Event Res", (digi->GetCharge() - fMCQA[address]) / fMCQA[address]);
342 }
343 if (Links == 2 && links[0].GetEntry() != links[1].GetEntry()) {
344 fQA->CreateHist("Inter Event Res", 400., -1., 1.);
345 fQA->Fill("Inter Event Res", (digi->GetCharge() - fMCQA[address]) / fMCQA[address]);
346 }
347 }
348 else {
349 fQA->CreateHist("E FN", 200, 0., 10.0);
350 fQA->Fill("E FN", digi->GetCharge());
351 fQA->CreateHist("E FN max", 512, -12., 500.);
352 fQA->Fill("E FN max", temp[fMaxBin]);
353 fQA->CreateHist("E FN MC max", 512, -12., 500., 200, 0., 50.);
354 fQA->Fill("E FN MC max", temp[fMaxBin], fMCQA[address]);
355 fQA->CreateHist("E FN res", 400., -1., 1.);
356 fQA->Fill("E FN res", digi->GetCharge() - fMCQA[address]);
357 fQA->CreateHist("E FN rel", 400., -1., 1.);
358 fQA->Fill("E FN rel", (digi->GetCharge() - fMCQA[address]) / fMCQA[address]);
359 fQA->CreateHist("Links Res FN", 400., -1., 1., 5, -0.5, 4.5);
360 fQA->Fill("Links Res FN", digi->GetCharge() - fMCQA[address], digiMatch->GetNofLinks());
361
362 fQA->CreateProfile("1 Link Prof FN", 32, 0., 32.);
363 for (size_t i = 0; i < fPulseBuffer[address].first.size(); i++) {
364 fQA->FillProfile("1 Link Prof FN", i, temp[i]);
365 }
366
367
368 if (fLinkQA[address].size() == 2 && fLinkQA[address][1]["Event"] != fLinkQA[address][0]["Event"]) {
369 fQA->CreateHist("Links QA time FN", 400., -1., 1., 500, 0., 2000.);
370 fQA->Fill("Links QA time FN", (digi->GetCharge() - fMCQA[address]) / fMCQA[address],
371 fLinkQA[address][1]["Time"] - fLinkQA[address][0]["Time"]);
372 fQA->CreateProfile("2 Link Prof FN", 32, 0., 32., 100, 0., 2000.);
373 for (size_t i = 0; i < fPulseBuffer[address].first.size(); i++) {
374 fQA->FillProfile("2 Link Prof FN", i, fLinkQA[address][1]["Time"] - fLinkQA[address][0]["Time"], temp[i]);
375 }
376 }
377 }
378 fShiftQA.erase(address);
379 fMCQA.erase(address);
380 fLinkQA.erase(address);
381 }
382
383
384 // digi->SetAddressModule(fModAddress); Not required anymore, now handled in the digi c'tor
385
386 if (trigger == 1) {
388 }
389 if (trigger == 0 && FNcall) {
391 }
392 if (trigger == 1 && MultiCall) {
394 }
395
396 //digi->SetMatch(digiMatch);
397 if (fDebug) {
398 fQA->CreateHist("Links", 10, -0.5, 9.5);
399 fQA->Fill("Links", digiMatch->GetNofLinks());
400 }
401 if (!FNcall && fPrintPulse)
402 std::cout << "main call charge: " << digi->GetCharge() << " col : " << col
403 << " lay : " << CbmTrdAddress::GetLayerId(address) << " trigger: " << trigger
404 << " time: " << digi->GetTime() << std::endl;
405 if (FNcall && fPrintPulse)
406 std::cout << "FN call charge: " << digi->GetCharge() << " col : " << col
407 << " lay : " << CbmTrdAddress::GetLayerId(address) << " trigger: " << trigger
408 << " time: " << digi->GetTime() << std::endl;
409
410 // if(!FNcall && MultiCall) std::cout<<"main call charge: "<<digi->GetCharge()<<" col : " << col<<" lay : " << CbmTrdAddress::GetLayerId(address)<<" trigger: " << trigger<<" time: " << digi->GetTime()<<std::endl;
411 // if(FNcall && MultiCall) std::cout<<"FN call charge: "<<digi->GetCharge()<<" col : " << col<<" lay : " << CbmTrdAddress::GetLayerId(address)<< " trigger: " << trigger<<" time: " << digi->GetTime()<<std::endl;
412
413 fDigiMap[address] = std::make_pair(digi, digiMatch);
414
415 delete fPulseBuffer[address].second;
416 fPulseBuffer.erase(address);
417
418 if (!FNcall && !MultiCall && trigger == 1) {
419 if (down) {
420 Int_t FNaddress = 0;
421 if (col >= 1)
423 CbmTrdAddress::GetModuleId(fModAddress), sec, row, col - 1);
424 Double_t timediff = TMath::Abs(fTimeBuffer[address] - fTimeBuffer[FNaddress]);
425 if (FNaddress != 0 && timediff < CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC))
426 ProcessPulseBuffer(FNaddress, true, false, true, false);
427 }
428
429 if (up) {
430 Int_t FNaddress = 0;
431 if (col < (ncols - 1))
433 sec, row, col + 1);
434 Double_t timediff = TMath::Abs(fTimeBuffer[address] - fTimeBuffer[FNaddress]);
435 if (FNaddress != 0 && timediff < CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC))
436 ProcessPulseBuffer(FNaddress, true, false, false, true);
437 }
438 }
439
440
441 if (!FNcall && MultiCall && trigger == 1) {
442 if (down) {
443 Int_t FNaddress = 0;
444 if (col >= 1)
446 CbmTrdAddress::GetModuleId(fModAddress), sec, row, col - 1);
447 // Double_t timediff = TMath::Abs(fTimeBuffer[address]-fTimeBuffer[FNaddress]);
448 if (FNaddress != 0) ProcessPulseBuffer(FNaddress, true, true, true, false);
449 }
450
451 if (up) {
452 Int_t FNaddress = 0;
453 if (col < (ncols - 1))
455 sec, row, col + 1);
456 // Double_t timediff = TMath::Abs(fTimeBuffer[address]-fTimeBuffer[FNaddress]);
457 if (FNaddress != 0) ProcessPulseBuffer(FNaddress, true, true, false, true);
458 }
459 }
460
461 fTimeBuffer.erase(address);
462}
463
464
465//_______________________________________________________________________________
466void CbmTrdModuleSimR::AddDigitoBuffer(Int_t address, Double_t charge, Double_t /*chargeTR*/, Double_t time,
467 Int_t trigger)
468{
469 Double_t weighting = charge;
471 TVector3 padPos, padPosErr;
472 fDigiPar->GetPadPosition(address, padPos, padPosErr);
473 Double_t distance = sqrt(pow(fXYZ[0] - padPos[0], 2) + pow(fXYZ[1] - padPos[1], 2));
474 weighting = 1. / distance;
475 }
476
477 //compare times of the buffer content with the actual time and process the buffer if collecttime is over
478 Bool_t eventtime = false;
479 if (time > 0.000) eventtime = true;
480 if (eventtime) CheckTime(address);
481
482 AddNoise(charge);
483
484 //fill digis in the buffer
485 CbmMatch* digiMatch = new CbmMatch();
486 digiMatch->AddLink(CbmLink(weighting, fPointId, fEventId, fInputId));
487
488 Int_t col = CbmTrdAddress::GetColumnId(address);
489 Int_t row = CbmTrdAddress::GetRowId(address);
490 Int_t sec = CbmTrdAddress::GetSectorId(address);
491 Int_t ncols = fDigiPar->GetNofColumns();
492 Int_t channel(0);
493 for (Int_t isec(0); isec < sec; isec++)
494 channel += ncols * fDigiPar->GetNofRowsInSector(isec);
495 channel += ncols * row + col;
496
497 auto triggertype = CbmTrdDigi::eTriggerType::kNTrg;
498 if (trigger == 1) triggertype = CbmTrdDigi::eTriggerType::kSelf;
499 if (trigger == 2) triggertype = CbmTrdDigi::eTriggerType::kNeighbor;
500 // std::cout<<charge*1e6<<" "<<fTimeBuffer[address]/CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC)<<std::endl;
501 CbmTrdDigi* digi =
502 new CbmTrdDigi(channel, fModAddress, charge * 1e6,
503 ULong64_t(time / CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC)), triggertype, 0);
504
505
506 //digi->SetMatch(digiMatch);
507 // printf("CbmTrdModuleSimR::AddDigitoBuffer(%10d)=%3d [%d] col[%3d] row[%d] sec[%d]\n", address, channel, fModAddress, col, row, sec);
508
509 fAnalogBuffer[address].push_back(std::make_pair(digi, digiMatch));
510 fTimeBuffer[address] = fCurrentTime;
511 if (!eventtime) ProcessBuffer(address);
512}
513
514//_______________________________________________________________________________
515void CbmTrdModuleSimR::AddDigitoPulseBuffer(Int_t address, Double_t reldrift, Double_t charge, Double_t /*chargeTR*/,
516 Double_t time, Int_t /*trigger*/, Int_t /*epoints*/, Int_t /*ipoint*/)
517{
518
519 Double_t weighting = charge * fepoints;
521 TVector3 padPos, padPosErr;
522 fDigiPar->GetPadPosition(address, padPos, padPosErr);
523 Double_t distance = sqrt(pow(fXYZ[0] - padPos[0], 2) + pow(fXYZ[1] - padPos[1], 2));
524 weighting = 1. / distance;
525 }
526
527 // if(!CbmTrdDigitizer::IsTimeBased() && fPointId-fLastPoint!=0) {CheckBuffer(true);CleanUp(true);}
528
529 if (CbmTrdDigitizer::IsTimeBased() && reldrift == 0.) {
530 Bool_t gotMulti = CheckMulti(address, fPulseBuffer[address].first);
531 fMultiBuffer.erase(address);
532 if (gotMulti) fMCBuffer[address].clear();
533 }
534
535 if (CbmTrdDigitizer::IsTimeBased() && reldrift == 0.) CheckTime(address);
536
537 fMultiBuffer[address] = std::make_pair(fMultiBuffer[address].first + charge, 0.);
538 if (fMCBuffer[address].size() > 0) {
539 fMCBuffer[address][0].push_back(fPointId);
540 fMCBuffer[address][1].push_back(fEventId);
541 fMCBuffer[address][2].push_back(fInputId);
542 }
543 else
544 for (auto ini = 0; ini < 3; ini++) {
545 std::vector<Int_t> v;
546 fMCBuffer[address].push_back(v);
547 }
548 std::vector<Double_t> pulse;
549 if (fTimeBuffer[address] > 0) {
550 pulse = fPulseBuffer[address].first;
551 if (pulse.size() < 32) return;
552 AddToPulse(address, charge, reldrift, pulse);
553 Bool_t found = false;
554 for (Int_t links = 0; links < fPulseBuffer[address].second->GetNofLinks(); links++)
555 if (fPulseBuffer[address].second->GetLink(links).GetIndex() == fPointId) found = true;
556 if (!found) {
557 fPulseBuffer[address].second->AddLink(CbmLink(weighting, fPointId, fEventId, fInputId));
558 if (fDebug) {
559 std::map<TString, Int_t> LinkQA;
560 LinkQA["Event"] = fEventId;
561 LinkQA["Point"] = fPointId;
562 LinkQA["Time"] = fEventTime;
563 LinkQA["Charge"] = charge * 1e6 * fepoints;
564 LinkQA["PosX"] = fQAPosition[0];
565 LinkQA["PosY"] = fQAPosition[1];
566 fLinkQA[address].push_back(LinkQA);
567 }
568 }
569 }
570 else {
571 fMCBuffer[address].clear();
572 pulse = MakePulse(charge, pulse, address);
573 fPulseBuffer[address].first = pulse;
574 CbmMatch* digiMatch = new CbmMatch();
575 digiMatch->AddLink(CbmLink(weighting, fPointId, fEventId, fInputId));
576
577 if (fDebug) {
578 std::vector<std::map<TString, Int_t>> vecLink;
579 std::map<TString, Int_t> LinkQA;
580 LinkQA["Event"] = fEventId;
581 LinkQA["Point"] = fPointId;
582 LinkQA["Time"] = fEventTime;
583 LinkQA["Charge"] = charge * 1e6 * fepoints;
584 LinkQA["PosX"] = fQAPosition[0];
585 LinkQA["PosY"] = fQAPosition[1];
586 vecLink.push_back(LinkQA);
587 fLinkQA[address] = vecLink;
588 }
589
590 fPulseBuffer[address].second = digiMatch;
591 fTimeBuffer[address] = time;
592 }
593
594 // if(!CbmTrdDigitizer::IsTimeBased() && finish) {CheckBuffer(true);CleanUp(true);}
595}
596
597
598std::vector<Double_t> CbmTrdModuleSimR::MakePulse(Double_t charge, std::vector<Double_t> pulse, Int_t address)
599{
600
601 Double_t sample[32];
602 for (Int_t i = 0; i < 32; i++)
603 sample[i] = 0;
604 // Double_t timeshift = gRandom->Uniform(0.,CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC));
605 Double_t timeshift =
607 if (fDebug) fQA->Fill("Shift", timeshift);
608 // Int_t shift=timeshift;
609 //fShiftQA[address]=shift;
610 fShiftQA[address] = timeshift;
611 for (Int_t i = 0; i < 32; i++) {
612 if (i < fPresamples) {
613 sample[i] = 0.;
614 continue;
615 }
616 if (fTimeShift)
617 sample[i] = fCalibration * charge * 1e6
619 if (!fTimeShift)
620 sample[i] = fCalibration * charge * 1e6
622 if (sample[i] > fClipLevel && fClipping) sample[i] = fClipLevel - 1; //clipping
623 }
624
625 for (Int_t i = 0; i < 32; i++) {
626 pulse.push_back(sample[i]);
627 }
628 AddCorrelatedNoise(pulse);
629
630 // if(fDebug && CheckTrigger(pulse) == 1) fMCQA[address]=charge*1e6;
631 if (fDebug) fMCQA[address] = charge * 1e6;
632 return pulse;
633}
634
635void CbmTrdModuleSimR::AddToPulse(Int_t address, Double_t charge, Double_t reldrift, std::vector<Double_t> pulse)
636{
637
638 Int_t comptrigger = CheckTrigger(pulse);
639 std::vector<Double_t> temppulse;
640 for (Int_t i = 0; i < 32; i++)
641 temppulse.push_back(pulse[i]);
642 Double_t dt = fCurrentTime - fTimeBuffer[address];
643 Int_t startbin = (dt + reldrift) / CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC);
644 Double_t timeshift =
645 ((Int_t)((dt + reldrift) * 10) % (Int_t)(CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC) * 10)) / 10.0;
646 if (startbin > 31 || dt < 0.) return;
647 if (fDebug) fMCQA[address] += charge * 1e6;
648
649 for (Int_t i = 0; i < 32; i++) {
650 if (i < fPresamples + startbin) {
651 pulse[i] = temppulse[i];
652 continue;
653 }
654 Int_t addtime = i - startbin - fPresamples;
655 pulse[i] += fCalibration * charge * 1e6
657 if (pulse[i] > fClipLevel && fClipping) pulse[i] = fClipLevel - 1; //clipping
658 }
659
660 // std::vector<Int_t> newpulse;
661 // for(Int_t i=0;i<32;i++){
662 // if(i < fPresamples) continue;
663 // if(fTimeShift){
664 // Int_t sample = fCalibration*charge*1e6*CalcResponse((i-fPresamples)*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC)+timeshift);
665 // if(sample > fClipLevel && fClipping) sample=fClipLevel-1; //clipping
666 // newpulse.push_back(sample);
667 // }
668 // if(!fTimeShift){
669 // Int_t sample = fCalibration*charge*1e6*CalcResponse((i-fPresamples)*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC));
670 // if(sample > fClipLevel && fClipping) sample=fClipLevel-1; //clipping
671 // newpulse.push_back(sample);
672 // }
673 // }
674
675 Int_t trigger = CheckTrigger(pulse);
676 if (comptrigger == 0 && trigger == 1) {
677 for (Int_t i = 0; i < 32; i++) {
678 if (i < fPresamples && startbin >= fPresamples) {
679 pulse[i] = temppulse[i + startbin - fPresamples];
680 continue;
681 }
682 if (i < fPresamples && startbin < fPresamples) {
683 pulse[i] = temppulse[i + startbin];
684 continue;
685 }
686 Int_t addtime = i - fPresamples;
687 Int_t shift = startbin + i;
688 if (fTimeShift) {
689 if (shift < 32)
690 pulse[i] = temppulse[shift]
691 + fCalibration * charge * 1e6
693 else
694 pulse[i] = fCalibration * charge * 1e6
696 if (pulse[i] > fClipLevel && fClipping) pulse[i] = fClipLevel - 1; //clipping
697 }
698 if (!fTimeShift) {
699 if (shift < 32)
700 pulse[i] = temppulse[shift]
701 + fCalibration * charge * 1e6
703 else
704 pulse[i] =
706 if (pulse[i] > fClipLevel && fClipping) pulse[i] = fClipLevel - 1; //clipping
707 }
708 }
709 fTimeBuffer[address] = fCurrentTime;
710 }
711
712 if (trigger == 2) {
713 fMultiBuffer[address].second = fCurrentTime;
714 }
715
716 fPulseBuffer[address].first = pulse;
717}
718
719Bool_t CbmTrdModuleSimR::CheckMulti(Int_t address, std::vector<Double_t> pulse)
720{
721
722 Bool_t processed = false;
723 Int_t trigger = CheckTrigger(pulse);
724 if (trigger == 2) {
725 processed = true;
726 Int_t maincol = CbmTrdAddress::GetColumnId(address);
727 Int_t row = CbmTrdAddress::GetRowId(address);
728 Int_t sec = CbmTrdAddress::GetSectorId(address);
729 Int_t shift = GetMultiBin(pulse);
730 Int_t ncols = fDigiPar->GetNofColumns();
731 // Double_t timeshift = gRandom->Uniform(0.,CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC));
732 Double_t timeshift =
734 / 10.0;
735
736 std::vector<Double_t> temppulse;
737 std::map<Int_t, std::vector<Double_t>> templow;
738 std::map<Int_t, std::vector<Double_t>> temphigh;
739
740 temppulse = pulse;
741 for (Int_t i = 0; i < 32; i++) {
742 if (fDebug) fQA->CreateHist("Multi self", 32, -0.5, 31.5, 512, -12., 500.);
743 if (fDebug) fQA->Fill("Multi self", i, pulse[i]);
744 if (i >= shift) pulse[i] = 0.;
745 }
746 fPulseBuffer[address].first = pulse;
747
748 // std::cout<<"multi time: " << fTimeBuffer[address]<<" col: "<<maincol<<std::endl;
749 // Double_t dt=TMath::Abs(fMultiBuffer[address].second - fTimeBuffer[address]);
750
751
752 Int_t col = maincol;
753 Int_t FNaddress = 0;
754 Double_t FNshift = 0;
755 std::vector<Double_t> FNpulse = fPulseBuffer[address].first;
756 if (col >= 1)
758 sec, row, col - 1);
759 Int_t FNtrigger = 1;
760
761 while (FNtrigger == 1 && FNaddress != 0) {
762 if (col >= 1)
764 CbmTrdAddress::GetModuleId(fModAddress), sec, row, col - 1);
765 else
766 break;
767 col--;
768 FNpulse = fPulseBuffer[FNaddress].first;
769 templow[FNaddress] = FNpulse;
770 FNshift =
772 for (Int_t i = 0; i < 32; i++) {
773 if (i >= FNshift) FNpulse[i] = 0.;
774 }
775 FNtrigger = CheckTrigger(FNpulse);
776 fPulseBuffer[FNaddress].first = FNpulse;
777 if (col == 0) break;
778
779 // std::cout<<"col: "<<col<<" "<< fTimeBuffer[FNaddress]<<" FNaddress: " << FNaddress<<" FNtrigger: "<< FNtrigger<<" samples: " << fPulseBuffer[FNaddress].first.size()<<" time: " << fTimeBuffer[FNaddress]<<std::endl;
780 }
781
782 col = maincol;
783 FNaddress = 0;
784 if (col < (ncols - 1))
786 sec, row, col + 1);
787 FNtrigger = 1;
788
789 while (FNtrigger == 1 && FNaddress != 0) {
790 if (col < (ncols - 1))
792 CbmTrdAddress::GetModuleId(fModAddress), sec, row, col + 1);
793 else
794 break;
795 col++;
796 FNpulse = fPulseBuffer[FNaddress].first;
797 temphigh[FNaddress] = FNpulse;
798 FNshift =
800 for (Int_t i = 0; i < 32; i++) {
801 if (i >= FNshift) FNpulse[i] = 0.;
802 }
803 FNtrigger = CheckTrigger(FNpulse);
804 fPulseBuffer[FNaddress].first = FNpulse;
805 if (col == ncols - 1) break;
806 }
807 ProcessPulseBuffer(address, false, true, true, true);
808
809 for (Int_t i = 0; i < 32; i++) {
810 if (i < fPresamples) {
811 pulse[i] = temppulse[shift + i - fPresamples];
812 continue;
813 }
814 if (shift + i - fPresamples < 32) {
815 if (fTimeShift)
816 pulse[i] =
817 temppulse[shift + i - fPresamples]
818 + fCalibration * fMultiBuffer[address].first * 1e6
820 if (!fTimeShift)
821 pulse[i] = temppulse[shift + i - fPresamples]
822 + fCalibration * fMultiBuffer[address].first * 1e6
824 }
825 else {
826 if (fTimeShift)
827 pulse[i] =
828 fCalibration * fMultiBuffer[address].first * 1e6
830 if (!fTimeShift)
831 pulse[i] = fCalibration * fMultiBuffer[address].first * 1e6
833 }
834
835 // if(shift+i < 32){
836 // if(fTimeShift) pulse[i]=temppulse[shift+i]+fCalibration*fMultiBuffer[address].first*1e6*CalcResponse(i*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC)+timeshift);
837 // if(!fTimeShift) pulse[i]=temppulse[shift+i]+fCalibration*fMultiBuffer[address].first*1e6*CalcResponse(i*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC));
838 // }
839 // else{
840 // if(fTimeShift) pulse[i]=fCalibration*fMultiBuffer[address].first*1e6*CalcResponse(i*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC)+timeshift);
841 // if(!fTimeShift) pulse[i]=fCalibration*fMultiBuffer[address].first*1e6*CalcResponse(i*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC));
842 // }
843 }
844 for (Int_t i = 0; i < 32; i++) {
845 if (fDebug) fQA->CreateHist("Multi self after", 32, -0.5, 31.5, 512, -12., 500.);
846 if (fDebug) fQA->Fill("Multi self after", i, pulse[i]);
847 }
848
849 fTimeBuffer[address] = fMultiBuffer[address].second;
850 fPulseBuffer[address].first = pulse;
851
852 if (col < ncols - 1)
854 sec, row, col + 1);
855 FNtrigger = 1;
856
857 while (FNtrigger == 1 && FNaddress != 0) {
858 if (col < (ncols - 1))
860 CbmTrdAddress::GetModuleId(fModAddress), sec, row, col + 1);
861 else
862 break;
863 col++;
864 FNpulse = fPulseBuffer[FNaddress].first;
865 for (Int_t i = 0; i < 32; i++) {
866 if (i < fPresamples && temphigh[FNaddress].size() > 0.) {
867 pulse[i] = temphigh[FNaddress][shift + i - fPresamples];
868 continue;
869 }
870 if (fTimeShift) {
871 if ((size_t) shift + i - fPresamples < temphigh[FNaddress].size())
872 FNpulse[i] =
873 temphigh[FNaddress][shift + i - fPresamples]
874 + fCalibration * fMultiBuffer[FNaddress].first * 1e6
876 else
877 FNpulse[i] =
878 fCalibration * fMultiBuffer[FNaddress].first * 1e6
880 if (FNpulse[i] > fClipLevel && fClipping) FNpulse[i] = fClipLevel - 1; //clipping
881 }
882 if (!fTimeShift) {
883 if ((size_t) shift + i - fPresamples < temphigh[FNaddress].size())
884 FNpulse[i] = temphigh[FNaddress][shift + i - fPresamples]
885 + fCalibration * fMultiBuffer[FNaddress].first * 1e6
887 else
888 FNpulse[i] = fCalibration * fMultiBuffer[FNaddress].first * 1e6
890 if (FNpulse[i] > fClipLevel && fClipping) FNpulse[i] = fClipLevel - 1; //clipping
891 }
892
893 // if(fTimeShift){
894 // if(shift+i<32 && temphigh[FNaddress].size()>0) FNpulse[i]=temphigh[FNaddress][shift+i]+fCalibration*fMultiBuffer[FNaddress].first*1e6*CalcResponse(i*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC)+timeshift);
895 // else FNpulse[i]=fCalibration*fMultiBuffer[FNaddress].first*1e6*CalcResponse(i*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC)+timeshift);
896 // if(FNpulse[i] > fClipLevel && fClipping) FNpulse[i]=fClipLevel-1; //clipping
897 // }
898 // if(!fTimeShift){
899 // if(shift+i<32 && temphigh[FNaddress].size()>0) FNpulse[i]=temphigh[FNaddress][shift+i]+fCalibration*fMultiBuffer[FNaddress].first*1e6*CalcResponse(i*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC));
900 // else FNpulse[i]=fCalibration*fMultiBuffer[FNaddress].first*1e6*CalcResponse(i*CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC));
901 // if(FNpulse[i] > fClipLevel && fClipping) FNpulse[i]=fClipLevel-1; //clipping
902 // }
903 }
904
905 fPulseBuffer[FNaddress].first = FNpulse;
906 CbmMatch* FNdigiMatch = new CbmMatch();
907 if (!fMCBuffer[FNaddress].empty())
908 for (UInt_t links = 0; links < fMCBuffer[FNaddress][0].size(); links++)
909 FNdigiMatch->AddLink(CbmLink(fMultiBuffer[FNaddress].first, fMCBuffer[FNaddress][0][links],
910 fMCBuffer[FNaddress][1][links], fMCBuffer[FNaddress][2][links]));
911 fPulseBuffer[FNaddress].second = FNdigiMatch;
912 if (fDebug) {
913 std::vector<std::map<TString, Int_t>> vecLink;
914 std::map<TString, Int_t> LinkQA;
915 LinkQA["Event"] = fEventId;
916 LinkQA["Point"] = fPointId;
917 LinkQA["Time"] = fEventTime;
918 LinkQA["Charge"] = fMultiBuffer[FNaddress].first * 1e6 * fepoints;
919 vecLink.push_back(LinkQA);
920 fLinkQA[FNaddress] = vecLink;
921 }
922
923 FNtrigger = CheckTrigger(FNpulse);
924
925 fTimeBuffer[FNaddress] = fMultiBuffer[address].second;
926 fMultiBuffer.erase(FNaddress);
927 if (col == ncols - 1) break;
928 }
929
930 if (col >= 1)
932 sec, row, col - 1);
933 FNtrigger = 1;
934
935 while (FNtrigger == 1 && FNaddress != 0) {
936 if (col >= 1)
938 CbmTrdAddress::GetModuleId(fModAddress), sec, row, col - 1);
939 else
940 break;
941 col--;
942 FNpulse = fPulseBuffer[FNaddress].first;
943 for (Int_t i = 0; i < 32; i++) {
944 if (i < fPresamples && templow[FNaddress].size() > 0.) {
945 pulse[i] = templow[FNaddress][shift + i - fPresamples];
946 continue;
947 }
948 if (fTimeShift) {
949 if ((size_t) shift + i - fPresamples < templow[FNaddress].size())
950 FNpulse[i] =
951 templow[FNaddress][shift + i - fPresamples]
952 + fCalibration * fMultiBuffer[FNaddress].first * 1e6
954 else
955 FNpulse[i] =
956 fCalibration * fMultiBuffer[FNaddress].first * 1e6
958 if (FNpulse[i] > fClipLevel && fClipping) FNpulse[i] = fClipLevel - 1; //clipping
959 }
960 if (!fTimeShift) {
961 if ((size_t) shift + i - fPresamples < templow[FNaddress].size())
962 FNpulse[i] = templow[FNaddress][shift + i - fPresamples]
963 + fCalibration * fMultiBuffer[FNaddress].first * 1e6
965 else
966 FNpulse[i] = fCalibration * fMultiBuffer[FNaddress].first * 1e6
968 if (FNpulse[i] > fClipLevel && fClipping) FNpulse[i] = fClipLevel - 1; //clipping
969 }
970 }
971 fPulseBuffer[FNaddress].first = FNpulse;
972 CbmMatch* FNdigiMatch = new CbmMatch();
973 if (!fMCBuffer[FNaddress].empty())
974 for (UInt_t links = 0; links < fMCBuffer[FNaddress][0].size(); links++)
975 FNdigiMatch->AddLink(CbmLink(fMultiBuffer[FNaddress].first, fMCBuffer[FNaddress][0][links],
976 fMCBuffer[FNaddress][1][links], fMCBuffer[FNaddress][2][links]));
977 fPulseBuffer[FNaddress].second = FNdigiMatch;
978 if (fDebug) {
979 std::vector<std::map<TString, Int_t>> vecLink;
980 std::map<TString, Int_t> LinkQA;
981 LinkQA["Event"] = fEventId;
982 LinkQA["Point"] = fPointId;
983 LinkQA["Time"] = fEventTime;
984 LinkQA["Charge"] = fMultiBuffer[FNaddress].first * 1e6 * fepoints;
985 vecLink.push_back(LinkQA);
986 fLinkQA[FNaddress] = vecLink;
987 }
988
989 FNtrigger = CheckTrigger(FNpulse);
990
991 fTimeBuffer[FNaddress] = fMultiBuffer[address].second;
992 fMultiBuffer.erase(FNaddress);
993 if (col == 0) break;
994 }
995
996 CbmMatch* digiMatch = new CbmMatch();
997 if (!fMCBuffer[address].empty())
998 for (UInt_t links = 0; links < fMCBuffer[address][0].size(); links++)
999 digiMatch->AddLink(CbmLink(fMultiBuffer[address].first, fMCBuffer[address][0][links],
1000 fMCBuffer[address][1][links], fMCBuffer[address][2][links]));
1001 fPulseBuffer[address].second = digiMatch;
1002 if (fDebug) {
1003 std::vector<std::map<TString, Int_t>> vecLink;
1004 std::map<TString, Int_t> LinkQA;
1005 LinkQA["Event"] = fEventId;
1006 LinkQA["Point"] = fPointId;
1007 LinkQA["Time"] = fEventTime;
1008 LinkQA["Charge"] = fMultiBuffer[address].first * 1e6 * fepoints;
1009 vecLink.push_back(LinkQA);
1010 fLinkQA[address] = vecLink;
1011 }
1012 }
1013
1014 return processed;
1015}
1016
1017Int_t CbmTrdModuleSimR::CheckTrigger(std::vector<Double_t> pulse)
1018{
1019
1020 Int_t slope = 0;
1021 Bool_t trigger = false;
1022 Bool_t falling = false;
1023 Bool_t multihit = false;
1024 for (size_t i = 0; i < pulse.size(); i++) {
1025 if (i < pulse.size() - 1) slope = pulse[i + 1] - pulse[i];
1026 if (slope >= fTriggerSlope && !trigger) trigger = true;
1027 if (slope < 0 && trigger) falling = true;
1028 if (slope >= fTriggerSlope && trigger && falling && (Int_t) i > fMaxBin) multihit = true;
1029 }
1030
1031 if (!trigger && !multihit) return 0;
1032 if (trigger && !multihit) return 1;
1033 if (trigger && multihit) return 2;
1034
1035 return 0;
1036}
1037
1038Int_t CbmTrdModuleSimR::GetMultiBin(std::vector<Double_t> pulse)
1039{
1040
1041 Int_t slope = 0;
1042 Int_t startbin = 0;
1043 Bool_t trigger = false;
1044 Bool_t falling = false;
1045 for (size_t i = 0; i < pulse.size(); i++) {
1046 if (i < 31) slope = pulse[i + 1] - pulse[i];
1047 if (slope >= fTriggerSlope && !trigger) trigger = true;
1048 if (slope < 0 && trigger) falling = true;
1049 if (slope >= fTriggerSlope && trigger && falling && (Int_t) i > fMaxBin) startbin = i;
1050 }
1051
1052 return startbin;
1053}
1054
1055
1056//_______________________________________________________________________________
1057Double_t CbmTrdModuleSimR::CalcPRF(Double_t x, Double_t W, Double_t h)
1058{
1059 Float_t K3 = 0.525;
1060 Double_t SqrtK3 = sqrt(K3);
1061
1062 return std::fabs(-1. / (2. * atan(SqrtK3))
1063 * (atan(SqrtK3 * tanh(TMath::Pi() * (-2. + SqrtK3) * (W + 2. * x) / (8. * h)))
1064 + atan(SqrtK3 * tanh(TMath::Pi() * (-2. + SqrtK3) * (W - 2. * x) / (8. * h)))));
1065}
1066
1067//_______________________________________________________________________________
1069{
1070
1071 if (fShapingOrder == 1)
1072 return (t / fTau) * TMath::Exp(-(t / fTau));
1073 else
1074 return (t / fTau) * (t / fTau) * TMath::Exp(-(t / fTau));
1075}
1076
1077//_______________________________________________________________________________
1078Bool_t CbmTrdModuleSimR::DistributeCharge(Double_t pointin[3], Double_t pointout[3], Double_t delta[3], Double_t pos[3],
1079 Int_t ipoints)
1080{
1081 if (fDistributionMode == 1) {
1082 for (Int_t i = 0; i < 3; i++) {
1083 pos[i] = pointin[i] + (0.01) * delta[i] + 0.95 * delta[i] / fepoints * ipoints;
1084 }
1085 }
1086
1087 if (fDistributionMode == 5) {
1088 for (Int_t i = 0; i < 3; i++) {
1089 pos[i] = pointin[i] + 0.5 * delta[i];
1090 }
1091 }
1092
1093
1094 //in development
1095 if (fDistributionMode == 2) {
1096 Double_t lastpos[3] = {pointin[0], pointin[1], pointin[2]};
1097 Double_t dist_gas = TMath::Sqrt(delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2]);
1098 if (ipoints > 0)
1099 for (Int_t i = 0; i < 3; i++)
1100 lastpos[i] = pos[i];
1101 Double_t roll = gRandom->Integer(100);
1102 Double_t s = (GetStep(dist_gas, roll) / dist_gas) / 2;
1103
1104 if (TMath::Abs(lastpos[0] + s * delta[0]) > fDigiPar->GetSizeX()
1105 || TMath::Abs(lastpos[1] + s * delta[1]) > fDigiPar->GetSizeY() || (lastpos[2] + s * delta[2]) >= pointout[2]) {
1106 for (Int_t i = 0; i < 3; i++) {
1107 pos[i] = lastpos[i];
1108 }
1109 return true;
1110 }
1111 for (Int_t i = 0; i < 3; i++)
1112 pos[i] = lastpos[i] + s * delta[i];
1113
1114 return true;
1115 }
1116
1117 if (fDistributionMode == 3) {
1118 for (Int_t i = 0; i < 3; i++) {
1119 pos[i] = pointin[i] + (0.5 + ipoints) * delta[i];
1120 }
1121 }
1122
1123
1124 return true;
1125}
1126
1127//_______________________________________________________________________________
1129{
1130 // if(!CbmTrdDigitizer::IsTimeBased()) fPulseSwitch = false;
1131 if (fPulseSwitch && !fMessageConverter->GetSetter()) {
1132 if (fDebug) fQA->CreateHist("Pulse", 32, -0.5, 31.5, 512, -12., 500.);
1133 if (fDebug) fQA->CreateHist("Shift", 63, -0.5, 62.5);
1134 if (fDebug) fQA->CreateHist("Multi Quote", 4, -0.5, 3.5);
1135 // if(fDebug) fMessageConverter->SetDebug(true);
1136 if (fDebug) fMessageConverter->SetQA(fQA);
1137 std::vector<Int_t> recomask;
1138 for (Int_t i = frecostart; i <= frecostop; i++)
1139 recomask.push_back(i);
1140 fMessageConverter->SetRecoMask(recomask);
1141 fMessageConverter->SetCalibration(fCalibration);
1142 fMessageConverter->SetShapingOrder(fShapingOrder);
1143 fMessageConverter->SetTau(fTau);
1144 fMessageConverter->SetPresamples(fPresamples);
1145 fMessageConverter->SetMaxBin(fMaxBin);
1146 fMessageConverter->SetMinBin(fMinBin);
1147 //fMessageConverter->SetLookup(1);
1148 TString dir = getenv("VMCWORKDIR");
1149 TString filename = dir + "/parameters/trd/FeatureExtractionLookup.root";
1150 // if(fDistributionMode == 1 && fepoints == 3) filename= dir + "/parameters/trd/Feature_mode1_3points.root";
1151 // if(fDistributionMode == 1 && fepoints == 5) filename= dir + "/parameters/trd/Feature_mode1_5points.root";
1152 // if(fDistributionMode == 4) filename= dir + "/parameters/trd/Feature_mode4.root";
1153 fMessageConverter->SetReadFile(filename.Data());
1154 fMessageConverter->Init();
1155 fMessageConverter->SetSetter(true);
1156 }
1157
1158
1159 // calculate current physical time
1160 fCurrentTime = time + point->GetTime(); //+ AddDrifttime(gRandom->Integer(240))*1000; //convert to ns;
1161 fEventTime = time;
1162
1163 const Double_t nClusterPerCm = 1.0;
1164 Double_t point_in[3] = {point->GetXIn(), point->GetYIn(), point->GetZIn()};
1165 Double_t point_out[3] = {point->GetXOut(), point->GetYOut(), point->GetZOut()};
1166 Double_t local_point_out[3]; // exit point coordinates in local module cs
1167 Double_t local_point_in[3]; // entrace point coordinates in local module cs
1168 gGeoManager->cd(GetPath());
1169 gGeoManager->MasterToLocal(point_in, local_point_in);
1170 gGeoManager->MasterToLocal(point_out, local_point_out);
1171 SetPositionMC(local_point_out);
1172
1173 for (Int_t i = 0; i < 3; i++)
1174 fQAPosition[i] = point_in[i];
1175 for (Int_t i = 0; i < 3; i++)
1176 fQAPos_out[i] = point_out[i];
1177
1178 //fCurrentTime -= 1e9 * (point_in[2] / 100)
1179 // / 2.99e8; // subtract time of flight to the detector layer
1180
1181 // General processing on the MC point
1182 Double_t ELoss(0.), ELossTR(0.), ELossdEdX(point->GetEnergyLoss());
1183 if (fRadiator && TR) {
1184 nofElectrons++;
1185 if (
1186 fRadiator->LatticeHit(
1187 point)) { // electron has passed lattice grid (or frame material) befor reaching the gas volume -> TR-photons have been absorbed by the lattice grid
1189 }
1190 else if (point_out[2] >= point_in[2]) { //electron has passed the radiator
1191 TVector3 mom;
1192 point->Momentum(mom);
1193 ELossTR = fRadiator->GetTR(mom);
1194 }
1195 }
1196 ELoss = ELossTR + ELossdEdX;
1197
1198 if (fDebug) fQA->CreateHist("E MC", 200, 0., 50.0);
1199 if (fDebug) fQA->Fill("E MC", ELoss * 1e6);
1200
1202
1203 Double_t cluster_pos[3]; // cluster position in local module coordinate system
1204 Double_t cluster_delta
1205 [3]; // vector pointing in MC-track direction with length of one path slice within chamber volume to creat n cluster
1206
1207 Double_t trackLength = 0;
1208
1209 for (Int_t i = 0; i < 3; i++) {
1210 cluster_delta[i] = (local_point_out[i] - local_point_in[i]);
1211 trackLength += cluster_delta[i] * cluster_delta[i];
1212 }
1213 trackLength = TMath::Sqrt(trackLength);
1214 Int_t nCluster =
1215 trackLength / nClusterPerCm + 0.9; // Track length threshold of minimum 0.1cm track length in gas volume
1216
1217 if (fnClusterConst > 0) {
1218 nCluster = fnClusterConst; // Set number of cluster to constant value
1219 }
1220
1221 if (nCluster < 1) {
1222 return kFALSE;
1223 }
1224 nCluster = 1;
1225
1226 for (Int_t i = 0; i < 3; i++) {
1227 cluster_delta[i] /= Double_t(nCluster);
1228 }
1229
1230 Double_t clusterELoss = ELoss / Double_t(nCluster);
1231 Double_t clusterELossTR = ELossTR / Double_t(nCluster);
1232
1233
1234 //to change the number of ionization points in the gas
1235 Int_t epoints = fepoints;
1236
1237 if (fDistributionMode == 3) {
1238 epoints = nCluster;
1239 }
1240 if (fDistributionMode == 5) {
1241 epoints = 1;
1242 }
1243
1244 //in development
1245 std::vector<Double_t> vec;
1246 std::pair<Int_t, std::vector<Double_t>> steps = std::make_pair(0, vec);
1247 if (fDistributionMode == 4) {
1248 Double_t dist_gas = TMath::Sqrt(cluster_delta[0] * cluster_delta[0] + cluster_delta[1] * cluster_delta[1]
1249 + cluster_delta[2] * cluster_delta[2]);
1250 steps = (GetTotalSteps(local_point_in, local_point_out, dist_gas));
1251 epoints = steps.first;
1252 fepoints = steps.first;
1253 if (fDebug) fQA->CreateHist("Ionization", 63, -0.5, 62.5);
1254 if (fDebug) fQA->Fill("Ionization", steps.first);
1255 if (fDebug) fQA->CreateHist("Dist Ionization", 63, -0.5, 62.5, 20, 0., 2.);
1256 if (fDebug) fQA->Fill("Dist Ionization", steps.first, dist_gas);
1257 }
1258
1259 if (epoints != 1) {
1260 clusterELoss = ELoss / epoints;
1261 clusterELossTR = ELossTR / epoints;
1262 }
1263
1264 if (!fPulseSwitch) {
1265 for (Int_t ipoints = 0; ipoints < epoints; ipoints++) {
1266 Bool_t dist = DistributeCharge(local_point_in, local_point_out, cluster_delta, cluster_pos, ipoints);
1267 if (!dist) break;
1268 if (fDistributionMode == 4)
1269 for (Int_t i = 0; i < 3; i++)
1270 cluster_pos[i] = steps.second[i + ipoints * 3];
1271
1272 if (fDigiPar->GetSizeX() < std::fabs(cluster_pos[0]) || fDigiPar->GetSizeY() < std::fabs(cluster_pos[1])) {
1273 printf("-> nC %i/%i x: %7.3f y: %7.3f \n", ipoints, nCluster - 1, cluster_pos[0], cluster_pos[1]);
1274 for (Int_t i = 0; i < 3; i++)
1275 printf(" (%i) | in: %7.3f + delta: %7.3f * cluster: %i/%i = "
1276 "cluster_pos: %7.3f out: %7.3f g_in:%f g_out:%f\n",
1277 i, local_point_in[i], cluster_delta[i], ipoints, (Int_t) nCluster, cluster_pos[i], local_point_out[i],
1278 point_in[i], point_out[i]);
1279 }
1280
1282 Int_t noiserate = gRandom->Uniform(0, 3); //still in development
1283 Double_t simtime = fCurrentTime;
1284 for (Int_t ndigi = 0; ndigi < noiserate; ndigi++) {
1285 NoiseTime(time);
1286 // ScanPadPlane(cluster_pos, gRandom->Gaus(0, fSigma_noise_keV * 1.E-6), 0,epoints,ipoints);
1287 }
1288 fCurrentTime = simtime;
1289 }
1290 ScanPadPlane(cluster_pos, 0., clusterELoss, clusterELossTR, epoints, ipoints);
1291 }
1292 }
1293
1294 Double_t driftcomp = 10000;
1295 Int_t start = -1;
1296 Double_t Ionizations[epoints][3];
1297 if (fPulseSwitch) {
1298 for (Int_t ipoints = 0; ipoints < epoints; ipoints++) {
1299 Bool_t dist = DistributeCharge(local_point_in, local_point_out, cluster_delta, cluster_pos, ipoints);
1300 if (!dist) break;
1301 if (fDistributionMode == 4)
1302 for (Int_t i = 0; i < 3; i++)
1303 cluster_pos[i] = steps.second[i + ipoints * 3];
1304 for (Int_t i = 0; i < 3; i++)
1305 Ionizations[ipoints][i] = cluster_pos[i];
1306
1307 if (fDigiPar->GetSizeX() < std::fabs(cluster_pos[0]) || fDigiPar->GetSizeY() < std::fabs(cluster_pos[1])) {
1308 printf("-> nC %i/%i x: %7.3f y: %7.3f \n", ipoints, nCluster - 1, cluster_pos[0], cluster_pos[1]);
1309 for (Int_t i = 0; i < 3; i++)
1310 printf(" (%i) | in: %7.3f + delta: %7.3f * cluster: %i/%i = "
1311 "cluster_pos: %7.3f out: %7.3f g_in:%f g_out:%f\n",
1312 i, local_point_in[i], cluster_delta[i], ipoints, (Int_t) nCluster, cluster_pos[i], local_point_out[i],
1313 point_in[i], point_out[i]);
1314 }
1315
1316 fDigiPar->ProjectPositionToNextAnodeWire(cluster_pos);
1317 Int_t relz = 239 - (cluster_pos[2] - local_point_in[2]) / 0.005;
1318 if (relz > 239 || relz < 0) relz = 239;
1319
1320 Double_t Drift_x =
1321 TMath::Abs(double(int(cluster_pos[0] * 1000000) % int(fDigiPar->GetAnodeWireSpacing() * 100)) / 100);
1322 // std::cout<<" pos: " << Drift_x<< " " << int(cluster_pos[0]*100000)<< " " << int(fDigiPar->GetAnodeWireSpacing()*100)<<" " << (int(cluster_pos[0]*100000) % int(fDigiPar->GetAnodeWireSpacing()*100))<<std::endl;
1323 if (TMath::Abs(AddDrifttime(Drift_x, cluster_pos[2])) < driftcomp
1324 && TMath::Abs(AddDrifttime(Drift_x, cluster_pos[2])) > 0.) {
1325 driftcomp = TMath::Abs(AddDrifttime(Drift_x, cluster_pos[2]));
1326 start = ipoints;
1327 }
1328 }
1329
1330 if (start == -1) return false;
1331 for (Int_t i = 0; i < 3; i++)
1332 cluster_pos[i] = Ionizations[start][i];
1333
1334 Int_t relz = 239 - (cluster_pos[2] - local_point_in[2]) / 0.005;
1335 if (relz > 239 || relz < 0) relz = 239;
1336 Double_t Drift_x =
1337 TMath::Abs(double(int(cluster_pos[0] * 1000000) % int(fDigiPar->GetAnodeWireSpacing() * 100)) / 100);
1338 Double_t reldrift = TMath::Abs(AddDrifttime(Drift_x, cluster_pos[2]) - driftcomp) * 1000;
1339 fDriftStart = driftcomp * 1000;
1341 if (reldrift < 250.) ScanPadPlane(cluster_pos, 0., clusterELoss, clusterELossTR, epoints, start);
1342 if (fPrintPulse) std::cout << std::endl;
1343
1344 for (Int_t ipoints = 0; ipoints < epoints; ipoints++) {
1345 if (ipoints == start) continue;
1346 for (Int_t i = 0; i < 3; i++)
1347 cluster_pos[i] = Ionizations[ipoints][i];
1348
1349 relz = 239 - (cluster_pos[2] - local_point_in[2]) / 0.005;
1350 if (relz > 239 || relz < 0) relz = 239;
1351 Drift_x = TMath::Abs(double(int(cluster_pos[0] * 1000000) % int(fDigiPar->GetAnodeWireSpacing() * 100)) / 100);
1352 reldrift = TMath::Abs(AddDrifttime(Drift_x, cluster_pos[2]) - driftcomp) * 1000;
1353 if (reldrift < 250.) ScanPadPlane(cluster_pos, reldrift, clusterELoss, clusterELossTR, epoints, ipoints);
1354 if (fPrintPulse) std::cout << std::endl;
1355 }
1356 }
1357
1358 fLastEventTime = time;
1362 return true;
1363}
1364
1365
1366//_______________________________________________________________________________
1367void CbmTrdModuleSimR::ScanPadPlane(const Double_t* local_point, Double_t reldrift, Double_t clusterELoss,
1368 Double_t clusterELossTR, Int_t epoints, Int_t ipoint)
1369{
1370 Int_t sectorId(-1), columnId(-1), rowId(-1);
1371 fDigiPar->GetPadInfo(local_point, sectorId, columnId, rowId);
1372 if (sectorId < 0 && columnId < 0 && rowId < 0) {
1373 return;
1374 }
1375 else {
1376 for (Int_t i = 0; i < sectorId; i++) {
1377 rowId += fDigiPar->GetNofRowsInSector(i); // local -> global row
1378 }
1379
1380 // for (Int_t i = 0; i < 3; i++) fQAPosition[i]=local_point[i];
1381
1382 Double_t displacement_x(0), displacement_y(0); //mm
1383 Double_t h = fDigiPar->GetAnodeWireToPadPlaneDistance();
1384 Double_t W(fDigiPar->GetPadSizeX(sectorId)), H(fDigiPar->GetPadSizeY(sectorId));
1385 fDigiPar->TransformToLocalPad(local_point, displacement_x, displacement_y);
1386
1387 Int_t maxRow(6);
1388 if (fnScanRowConst > 0) maxRow = fnScanRowConst;
1389
1390 Int_t startRow(rowId - maxRow / 2);
1391 Int_t secRow(-1), targCol(-1), targRow(-1), targSec(-1), address(-1), fnRow(fDigiPar->GetNofRows()),
1392 fnCol(fDigiPar->GetNofColumns());
1393
1394 for (Int_t iRow = startRow; iRow <= rowId + maxRow / 2; iRow++) {
1395 Int_t iCol = columnId;
1396 if (((iCol >= 0) && (iCol <= fnCol - 1)) && ((iRow >= 0) && (iRow <= fnRow - 1))) { // real adress
1397 targSec = fDigiPar->GetSector(iRow, secRow);
1399 }
1400 else {
1401 targRow = iRow;
1402 targCol = iCol;
1403 if (iCol < 0) {
1404 targCol = 0;
1405 }
1406 else if (iCol > fnCol - 1) {
1407 targCol = fnCol - 1;
1408 }
1409 if (iRow < 0) {
1410 targRow = 0;
1411 }
1412 else if (iRow > fnRow - 1) {
1413 targRow = fnRow - 1;
1414 }
1415
1416 targSec = fDigiPar->GetSector(targRow, secRow);
1417 address =
1419 }
1420
1421 Bool_t print = false;
1422
1423 //distribute the mc charge fraction over the channels wit the PRF
1424 Float_t chargeFraction = 0;
1425 Float_t ch = 0;
1426 Float_t tr = 0;
1427
1428 // std::cout<<" prf half: " << CalcPRF(0 * W , W, h)<<std::endl;
1429 // std::cout<<" prf half -1 : " << CalcPRF(-1 * W , W, h)<<std::endl;
1430 // std::cout<<" prf half +1: " << CalcPRF(1 * W , W, h)<<std::endl;
1431 chargeFraction =
1432 CalcPRF((iCol - columnId) * W - displacement_x, W, h) * CalcPRF((iRow - rowId) * H - displacement_y, H, h);
1433
1434 ch = chargeFraction * clusterELoss;
1435 tr = chargeFraction * clusterELossTR;
1436
1437 Bool_t lowerend = false;
1438 Bool_t upperend = false;
1439 Int_t collow = 1;
1440 Int_t colhigh = 1;
1441
1442 if (fDebug) fQA->CreateHist("E self MC", 200, 0., 50.0);
1443 if (fDebug) fQA->CreateHist("E FN MC", 200, 0., 10.0);
1444
1445 if (ch >= (fMinimumChargeTH / epoints)) {
1446 if (!CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && !print)
1447 AddDigi(address, ch, tr, fCurrentTime, Int_t(1));
1448 if (!CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && print) {
1449 AddDigi(address, ch, tr, fCurrentTime, Int_t(1));
1450 std::cout << " time: " << fCurrentTime << " col: " << iCol << " row: " << iRow - rowId
1451 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 1 " << std::endl;
1452 }
1454 if (fPulseSwitch) {
1455 if (fDebug) fQA->Fill("E self MC", ch * epoints * 1e6);
1456 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(1), epoints, ipoint);
1457 }
1458 if (fPulseSwitch && print) {
1459 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(1), epoints, ipoint);
1460 std::cout << " time: " << fCurrentTime << " col: " << iCol << " row: " << iRow - rowId
1461 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 1 " << std::endl;
1462 }
1463
1464
1465 while (!lowerend) {
1466 if ((((iCol - collow) >= 0) && ((iCol - collow) <= fnCol - 1))
1467 && ((iRow >= 0) && (iRow <= fnRow - 1))) { // real adress
1468 targSec = fDigiPar->GetSector(iRow, secRow);
1470 iCol - collow);
1471 }
1472 else {
1473 break;
1474 }
1475
1476 chargeFraction = CalcPRF(((iCol - collow) - columnId) * W - displacement_x, W, h)
1477 * CalcPRF((iRow - rowId) * H - displacement_y, H, h);
1478 ch = chargeFraction * clusterELoss;
1479 tr = chargeFraction * clusterELossTR;
1480
1481
1482 if (ch >= (fMinimumChargeTH / epoints) && !CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && !print) {
1483 AddDigi(address, ch, tr, fCurrentTime, Int_t(1));
1484 collow++;
1485 }
1486 if (ch < (fMinimumChargeTH / epoints) && !CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && !print) {
1487 AddDigi(address, ch, tr, fCurrentTime, Int_t(2));
1488 lowerend = true;
1489 }
1490 if (ch >= (fMinimumChargeTH / epoints) && !CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && print) {
1491 std::cout << " time: " << fCurrentTime << " col: " << iCol - collow << " row: " << iRow - rowId
1492 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 1 " << std::endl;
1493 AddDigi(address, ch, tr, fCurrentTime, Int_t(1));
1494 collow++;
1495 }
1496 if (ch < (fMinimumChargeTH / epoints) && !CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && print) {
1497 std::cout << " time: " << fCurrentTime << " col: " << iCol - collow << " row: " << iRow - rowId
1498 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 0 " << std::endl;
1499 AddDigi(address, ch, tr, fCurrentTime, Int_t(2));
1500 lowerend = true;
1501 }
1502 if (ch >= (fMinimumChargeTH / epoints) && CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch) {
1503 AddDigitoBuffer(address, ch, tr, fCurrentTime, Int_t(1));
1504 collow++;
1505 }
1506 if (ch < (fMinimumChargeTH / epoints) && CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch) {
1507 AddDigitoBuffer(address, ch, tr, fCurrentTime, Int_t(2));
1508 lowerend = true;
1509 }
1510 if (ch >= (fMinimumChargeTH / epoints) && fPulseSwitch && !print) {
1511 if (fDebug) fQA->Fill("E self MC", ch * epoints * 1e6);
1512 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(1), epoints, ipoint);
1513 collow++;
1514 }
1515 if (ch < (fMinimumChargeTH / epoints) && fPulseSwitch && !print) {
1516 if (fDebug) fQA->Fill("E FN MC", ch * epoints * 1e6);
1517 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(2), epoints, ipoint);
1518 lowerend = true;
1519 }
1520
1521 if (ch >= (fMinimumChargeTH / epoints) && fPulseSwitch && print) {
1522 std::cout << " time: " << fCurrentTime << " col: " << iCol - collow << " row: " << iRow - rowId
1523 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 1 " << std::endl;
1524 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(1), epoints, ipoint);
1525 collow++;
1526 }
1527 if (ch < (fMinimumChargeTH / epoints) && fPulseSwitch && print) {
1528 std::cout << " time: " << fCurrentTime << " col: " << iCol - collow << " row: " << iRow - rowId
1529 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 0 " << std::endl;
1530 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(2), epoints, ipoint);
1531 lowerend = true;
1532 }
1533 }
1534
1535 while (!upperend) {
1536
1537 if ((((iCol + colhigh) >= 0) && ((iCol + colhigh) <= fnCol - 1))
1538 && ((iRow >= 0) && (iRow <= fnRow - 1))) { // real adress
1539 targSec = fDigiPar->GetSector(iRow, secRow);
1541 iCol + colhigh);
1542 }
1543 else {
1544 break;
1545 }
1546
1547
1548 chargeFraction = CalcPRF(((iCol + colhigh) - columnId) * W - displacement_x, W, h)
1549 * CalcPRF((iRow - rowId) * H - displacement_y, H, h);
1550 ch = chargeFraction * clusterELoss;
1551 tr = chargeFraction * clusterELossTR;
1552
1553 if (ch >= (fMinimumChargeTH / epoints) && !CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && !print) {
1554 AddDigi(address, ch, tr, fCurrentTime, Int_t(1));
1555 colhigh++;
1556 }
1557 if (ch < (fMinimumChargeTH / epoints) && !CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && !print) {
1558 AddDigi(address, ch, tr, fCurrentTime, Int_t(2));
1559 upperend = true;
1560 }
1561 if (ch >= (fMinimumChargeTH / epoints) && !CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && print) {
1562 std::cout << " time: " << fCurrentTime << " col: " << iCol + colhigh << " row: " << iRow - rowId
1563 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 1 " << std::endl;
1564 AddDigi(address, ch, tr, fCurrentTime, Int_t(1));
1565 colhigh++;
1566 }
1567 if (ch < (fMinimumChargeTH / epoints) && !CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch && print) {
1568 std::cout << " time: " << fCurrentTime << " col: " << iCol + colhigh << " row: " << iRow - rowId
1569 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 0 " << std::endl;
1570 AddDigi(address, ch, tr, fCurrentTime, Int_t(2));
1571 upperend = true;
1572 }
1573 if (ch >= (fMinimumChargeTH / epoints) && CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch) {
1574 AddDigitoBuffer(address, ch, tr, fCurrentTime, Int_t(1));
1575 colhigh++;
1576 }
1577 if (ch < (fMinimumChargeTH / epoints) && CbmTrdDigitizer::IsTimeBased() && !fPulseSwitch) {
1578 AddDigitoBuffer(address, ch, tr, fCurrentTime, Int_t(2));
1579 upperend = true;
1580 }
1581 if (ch >= (fMinimumChargeTH / epoints) && fPulseSwitch && !print) {
1582 if (fDebug) fQA->Fill("E self MC", ch * epoints * 1e6);
1583 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(1), epoints, ipoint);
1584 colhigh++;
1585 }
1586 if (ch < (fMinimumChargeTH / epoints) && fPulseSwitch && !print) {
1587 if (fDebug) fQA->Fill("E FN MC", ch * epoints * 1e6);
1588 if (ipoint == epoints - 1 && epoints > 1)
1589 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(2), epoints, ipoint);
1590 if (ipoint != epoints - 1 && epoints > 1)
1591 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(2), epoints, ipoint);
1592 if (epoints == 1) AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(2), epoints, ipoint);
1593 upperend = true;
1594 }
1595
1596 if (ch >= (fMinimumChargeTH / epoints) && fPulseSwitch && print) {
1597 std::cout << " time: " << fCurrentTime << " col: " << iCol + colhigh << " row: " << iRow - rowId
1598 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 1 " << std::endl;
1599 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(1), epoints, ipoint);
1600 colhigh++;
1601 }
1602 if (ch < (fMinimumChargeTH / epoints) && fPulseSwitch && print) {
1603 std::cout << " time: " << fCurrentTime << " col: " << iCol + colhigh << " row: " << iRow - rowId
1604 << " secrow: " << secRow << " charge: " << ch * 1e6 << " 0 " << std::endl;
1605 if (ipoint == epoints - 1 && epoints > 1)
1606 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(2), epoints, ipoint);
1607 if (ipoint != epoints - 1 && epoints > 1)
1608 AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(2), epoints, ipoint);
1609 if (epoints == 1) AddDigitoPulseBuffer(address, reldrift, ch, tr, fCurrentTime, Int_t(2), epoints, ipoint);
1610 upperend = true;
1611 }
1612 }
1613
1614 if (print) std::cout << std::endl;
1615
1616 } //if charge > trigger
1617 } //for rows
1618 }
1619}
1620
1621
1622//_______________________________________________________________________________
1624{
1627
1628 if (!fDigiPar) {
1629 LOG(warn) << GetName() << "::SetAsicPar : No Digi params for module " << fModAddress
1630 << ". Try calling first CbmTrdModSim::SetDigiPar.";
1631 return;
1632 }
1633
1634 if (fAsicPar) {
1635 LOG(warn) << GetName() << "::SetAsicPar : The list for module " << fModAddress << " already initialized.";
1636 return;
1637 }
1638 fAsicPar = new CbmTrdParModAsic();
1639 CbmTrdParSpadic* asic(NULL);
1640
1641 Int_t iFebGroup = 0; // 1; 2; // normal, super, ultimate
1642 Int_t gRow[3] = {2, 2, 2}; // re-ordering on the feb -> same mapping for normal and super
1643 Int_t gCol[3] = {16, 8, 4}; // re-ordering on the feb -> same mapping for normal and super
1644 Double_t xAsic = 0; // x position of Asic
1645 Double_t yAsic = 0; // y position of Asic
1646
1647 Int_t rowId(0), isecId(0), irowId(0), iAsic(0);
1648 for (Int_t s = 0; s < fDigiPar->GetNofSectors(); s++) {
1649 for (Int_t r = 0; r < fDigiPar->GetNofRowsInSector(s); r++) {
1650 for (Int_t c = 0; c < fDigiPar->GetNofColumnsInSector(s); c++) {
1651 // ultimate density 6 rows, 5 pads
1652 // super density 4 rows, 8 pads
1653 // normal density 2 rows, 16 pads
1654 if ((rowId % gRow[iFebGroup]) == 0) {
1655 if ((c % gCol[iFebGroup]) == 0) {
1656 xAsic = c + gCol[iFebGroup] / 2.;
1657 yAsic = r + gRow[iFebGroup] / 2.;
1658
1659 Double_t local_point[3];
1660 Double_t padsizex = fDigiPar->GetPadSizeX(s);
1661 Double_t padsizey = fDigiPar->GetPadSizeY(s);
1662
1663 // calculate position in sector coordinate system
1664 // with the origin in the lower left corner (looking upstream)
1665 local_point[0] = ((Int_t)(xAsic + 0.5) * padsizex);
1666 local_point[1] = ((Int_t)(yAsic + 0.5) * padsizey);
1667
1668 // calculate position in module coordinate system
1669 // with the origin in the lower left corner (looking upstream)
1670 local_point[0] += fDigiPar->GetSectorBeginX(s);
1671 local_point[1] += fDigiPar->GetSectorBeginY(s);
1672
1673 // local_point[i] must be >= 0 at this point Double_t local_point[3];
1674 Double_t fDx(GetDx()), fDy(GetDy());
1675 asic = new CbmTrdParSpadic(iAsic, iFebGroup, local_point[0] - fDx, local_point[1] - fDy);
1676 fAsicPar->SetAsicPar(asic);
1677 if (local_point[0] > 2 * fDx) {
1678 LOG(error) << "CbmTrdModuleSimR::SetAsicPar: asic position x=" << local_point[0]
1679 << " is out of bounds [0," << 2 * fDx << "]!";
1680 fDigiPar->Print("all");
1681 }
1682 if (local_point[1] > 2 * fDy) {
1683 LOG(error) << "CbmTrdModuleSimR::SetAsicPar: asic position y=" << local_point[1]
1684 << " is out of bounds [0," << 2 * fDy << "]!";
1685 fDigiPar->Print("all");
1686 }
1687 for (Int_t ir = rowId; ir < rowId + gRow[iFebGroup]; ir++) {
1688 for (Int_t ic = c; ic < c + gCol[iFebGroup]; ic++) {
1689 if (ir >= fDigiPar->GetNofRows())
1690 LOG(error) << "CbmTrdModuleSimR::SetAsicPar: ir " << ir << " is out of bounds!";
1691 if (ic >= fDigiPar->GetNofColumns())
1692 LOG(error) << "CbmTrdModuleSimR::SetAsicPar: ic " << ic << " is out of bounds!";
1693 isecId = fDigiPar->GetSector((Int_t) ir, irowId);
1694 asic->SetChannelAddress(
1696 //s, ir, ic));//new
1697 if (false)
1698 printf(" M:%10i(%4i) s: %i irowId: %4i ic: "
1699 "%4i r: %4i c: %4i address:%10i\n",
1700 fModAddress, CbmTrdAddress::GetModuleId(fModAddress), isecId, irowId, ic, r, c,
1701 CbmTrdAddress::GetAddress(fLayerId, fModAddress, isecId, irowId, ic));
1702 }
1703 }
1704 iAsic++; // next Asic
1705 }
1706 }
1707 }
1708 rowId++;
1709 }
1710 }
1711
1712 // Self Test
1713 for (Int_t s = 0; s < fDigiPar->GetNofSectors(); s++) {
1714 const Int_t nRow = fDigiPar->GetNofRowsInSector(s);
1715 const Int_t nCol = fDigiPar->GetNofColumnsInSector(s);
1716 for (Int_t r = 0; r < nRow; r++) {
1717 for (Int_t c = 0; c < nCol; c++) {
1719 if (fAsicPar->GetAsicAddress(channelAddress) == -1)
1720 LOG(error) << "CbmTrdModuleSimR::SetAsicPar: Channel address:" << channelAddress
1721 << " is not or multible initialized in module " << fModAddress
1722 << "(ID:" << CbmTrdAddress::GetModuleId(fModAddress) << ")"
1723 << "(s:" << s << ", r:" << r << ", c:" << c << ")";
1724 }
1725 }
1726 }
1727}
1728
1729
1730//_______________________________________________________________________________
1731void CbmTrdModuleSimR::SetNoiseLevel(Double_t sigma_keV) { fSigma_noise_keV = sigma_keV; }
1732
1733//_______________________________________________________________________________
1735
1736//_______________________________________________________________________________
1737void CbmTrdModuleSimR::SetSpadicResponse(Double_t calibration, Double_t tau)
1738{
1739
1740 fCalibration = calibration;
1741 fTau = tau;
1742 Double_t sum = 0;
1743 for (auto i = frecostart; i <= frecostop; i++)
1745 fEReco = sum;
1746}
1747
1748//_______________________________________________________________________________
1750{
1751
1752 if (mode == 1) {
1753 frecostart = 2;
1754 frecostop = 8;
1755 }
1756 if (mode == 2) {
1757 frecostart = 2;
1758 frecostop = 6;
1759 }
1760}
1761
1762
1763//_______________________________________________________________________________
1764void CbmTrdModuleSimR::SetPulseMode(Bool_t pulsed = true) { fPulseSwitch = pulsed; }
1765
1766
1767//_______________________________________________________________________________
1769{
1770 if (row % 2 == 0) row += 1;
1771 fnScanRowConst = row;
1772}
1773
1774
1775//_______________________________________________________________________________
1776Double_t CbmTrdModuleSimR::AddNoise(Double_t charge)
1777{
1778
1780 Double_t noise = gRandom->Gaus(0,
1781 fSigma_noise_keV * 1e-6); // keV->GeV // add only once per digi and event noise !!!
1782 charge += noise; // resulting charge can be < 0 -> possible problems with position reconstruction
1783 return charge;
1784 }
1785 else
1786 return 0.;
1787}
1788
1789//_______________________________________________________________________________
1791{
1792
1793
1795 if (fAdcNoise == 0) return 0.;
1796 Int_t noise = gRandom->Gaus(0, fAdcNoise);
1797 return noise;
1798 // return 0;
1799 }
1800 else
1801 return 0.;
1802}
1803
1804std::vector<Double_t> CbmTrdModuleSimR::AddCorrelatedNoise(std::vector<Double_t> pulse)
1805{
1806
1807 //dummy for now
1808 return pulse;
1809
1810 // for(size_t i=0;i<pulse.size();i++){
1811 // pulse[i] += TMath::Sin(i)*5;
1812 // }
1813
1814 // return pulse;
1815}
1816
1817
1818//_______________________________________________________________________________
1819Int_t CbmTrdModuleSimR::AddCrosstalk(Double_t address, Int_t i, Int_t sec, Int_t row, Int_t col, Int_t ncols)
1820{
1821
1822 Double_t cross = 0.;
1823 if (fAddCrosstalk && fPulseSwitch) {
1824 Int_t FNaddress = 0;
1825 if (col >= 1)
1827 sec, row, col - 1);
1828 if (fPulseBuffer[FNaddress].first.size() >= (size_t) i + 1)
1829 cross += fPulseBuffer[FNaddress].first[i] * fCrosstalkLevel;
1830
1831 FNaddress = 0;
1832 if (col < (ncols - 1))
1834 sec, row, col + 1);
1835 if (fPulseBuffer[FNaddress].first.size() >= (size_t) i + 1)
1836 cross += fPulseBuffer[FNaddress].first[i] * fCrosstalkLevel;
1837 }
1838
1839 return cross;
1840}
1841
1842//_______________________________________________________________________________
1844{
1845
1846
1847 std::map<Int_t, Double_t>::iterator timeit;
1848 std::vector<Int_t> toBeErased;
1849
1850 Bool_t done = false;
1851
1852 while (!done) {
1853 done = true;
1854 for (timeit = fTimeBuffer.begin(); timeit != fTimeBuffer.end(); timeit++) {
1855 Int_t add = timeit->first;
1856 if (fCurrentTime < fTimeBuffer[add]) continue;
1857 Double_t dt = fCurrentTime - fTimeBuffer[add];
1858 if ((dt < fCollectTime || dt == fCurrentTime) && !EB) continue;
1859 // if(!fPulseSwitch) {ProcessBuffer(add);fTimeBuffer.erase(add);}
1860 if (!fPulseSwitch) {
1861 ProcessBuffer(add);
1862 toBeErased.push_back(add);
1863 }
1864 if (fPulseSwitch) {
1865 std::vector<Double_t> pulse;
1866 pulse = fPulseBuffer[add].first;
1867
1868 if (CheckTrigger(pulse) == 1 && EB) {
1869 ProcessPulseBuffer(add, false, false, true, true);
1870 break;
1871 }
1872
1873
1874 if (CheckTrigger(pulse) == 1 && !EB) {
1875 ProcessPulseBuffer(add, false, false, true, true);
1876 done = false;
1877 break;
1878 }
1879
1880 if (fPrintPulse) std::cout << std::endl;
1881 }
1882 }
1883 }
1884
1885 for (auto& address : toBeErased) {
1886 fTimeBuffer.erase(address);
1887 }
1888}
1889
1890//_______________________________________________________________________________
1892{
1893 Bool_t closeTS(kFALSE);
1894 if (fTimeSlice && (time - fTimeSlice->GetEndTime()) > -1000.) closeTS = true;
1895
1896 // if(!CbmTrdDigitizer::IsTimeBased()) return 0;
1897 //process channels before timeslice ending and release memory
1898 // if(closeTS && CbmTrdDigitizer::IsTimeBased()){
1899 if (closeTS || time == 0) {
1900 std::map<Int_t, Double_t>::iterator timeit;
1901 Bool_t done = false;
1902
1903 while (!done) {
1904 done = true;
1905 for (timeit = fTimeBuffer.begin(); timeit != fTimeBuffer.end(); timeit++) {
1906 Int_t add = timeit->first;
1907 if (!fPulseSwitch) {
1908 ProcessBuffer(add);
1909 }
1910 if (fPulseSwitch) {
1911 std::vector<Double_t> pulse;
1912 pulse = fPulseBuffer[add].first;
1913
1914 if (CheckTrigger(pulse) == 1) {
1915 ProcessPulseBuffer(add, false, false, true, true);
1916 done = false;
1917 break;
1918 }
1919
1920 if (fPrintPulse) std::cout << std::endl;
1921 }
1922 }
1923 }
1924 std::map<Int_t, std::pair<std::vector<Double_t>, CbmMatch*>>::iterator itBuffer;
1925 for (itBuffer = fPulseBuffer.begin(); itBuffer != fPulseBuffer.end(); itBuffer++) {
1926 if (fPulseBuffer[itBuffer->first].second) delete fPulseBuffer[itBuffer->first].second;
1927 }
1928 fPulseBuffer.clear();
1929 fTimeBuffer.clear();
1930 fMultiBuffer.clear();
1931 fMCBuffer.clear();
1932 }
1933 return 0;
1934}
1935
1936
1937//_______________________________________________________________________________
1939{
1940
1941 std::map<Int_t, Double_t>::iterator timeit;
1942 // clean up
1943 std::vector<Int_t> erase_list;
1944
1945 if (fPulseSwitch) {
1946 for (timeit = fTimeBuffer.begin(); timeit != fTimeBuffer.end(); timeit++) {
1947 Int_t add = timeit->first;
1948 Double_t dt = fCurrentTime - fTimeBuffer[add];
1949 if (fTimeSlice) {
1950 if (fTimeBuffer[add] < fTimeSlice->GetStartTime()) {
1951 erase_list.push_back(add);
1952 continue;
1953 }
1954 }
1955 if ((dt < fCollectTime || dt == fCurrentTime) && !EB) continue;
1956 erase_list.push_back(add);
1957 }
1958 }
1959
1960 // //release memory for all non-triggered channels after signal collection time
1961 for (UInt_t i = 0; i < erase_list.size(); i++) {
1962 if (fPulseBuffer[erase_list[i]].second) delete fPulseBuffer[erase_list[i]].second;
1963 fPulseBuffer.erase(erase_list[i]);
1964 fTimeBuffer.erase(erase_list[i]);
1965 fMultiBuffer.erase(erase_list[i]);
1966 fMCBuffer.erase(erase_list[i]);
1967 }
1968}
1969
1970//_______________________________________________________________________________
1972{
1973
1974 //compare last entry in the actual channel with the current time
1975 std::map<Int_t, Double_t>::iterator timeit;
1976 Double_t dt = fCurrentTime - fTimeBuffer[address];
1977 // std::cout<<" dt: " << dt<<std::endl;
1978 Bool_t go = false;
1979 if (fCurrentTime > fTimeBuffer[address] && dt > 0.0000000) {
1980 if (dt > fCollectTime && dt != fCurrentTime && !fPulseSwitch) {
1981 ProcessBuffer(address);
1982 fTimeBuffer.erase(address);
1983 go = true;
1984 }
1985 // if(dt>fCollectTime && dt!=fCurrentTime && fPulseSwitch) {ProcessPulseBuffer(address,false,false);std::cout<<" ------ " <<std::endl;go=true;}
1986 if (dt > fCollectTime && dt != fCurrentTime && fPulseSwitch) {
1987 //ProcessPulseBuffer(address,false,false,true,true);
1988 go = true;
1989 if (fPrintPulse) std::cout << std::endl;
1990 }
1991 }
1992
1993 if (go && fPulseSwitch) {
1994 CheckBuffer(false);
1995 CleanUp(false);
1996 }
1997 if (go && !fPulseSwitch) CheckBuffer();
1998}
1999
2000//_______________________________________________________________________________
2001void CbmTrdModuleSimR::NoiseTime(ULong64_t eventTime) { fCurrentTime = gRandom->Uniform(fLastEventTime, eventTime); }
2002
2003//_______________________________________________________________________________
2004Double_t CbmTrdModuleSimR::AddDrifttime(Double_t x, Double_t z)
2005{
2006
2007 if (fDebug) fQA->CreateHist("All drift", 300, 0., 300.);
2008 if (fDebug) fQA->Fill("All drift", fDriftTime->GetBinContent(fDriftTime->FindBin(x, z)) * 1000);
2009
2010 return fDriftTime->GetBinContent(fDriftTime->FindBin(x, z));
2011}
2012
2013//_______________________________________________________________________________
2015{
2016 Double_t drifttime[241] = {
2017 0.11829, 0.11689, 0.11549, 0.11409, 0.11268, 0.11128, 0.10988, 0.10847, 0.10707, 0.10567, 0.10427,
2018 0.10287, 0.10146, 0.10006, 0.09866, 0.09726, 0.095859, 0.094459, 0.09306, 0.091661, 0.090262, 0.088865,
2019 0.087467, 0.086072, 0.084677, 0.083283, 0.08189, 0.080499, 0.07911, 0.077722, 0.076337, 0.074954, 0.073574,
2020 0.072197, 0.070824, 0.069455, 0.06809, 0.066731, 0.065379, 0.064035, 0.0627, 0.061376, 0.060063, 0.058764,
2021 0.05748, 0.056214, 0.054967, 0.053743, 0.052544, 0.051374, 0.05024, 0.049149, 0.048106, 0.047119, 0.046195,
2022 0.045345, 0.044583, 0.043925, 0.043403, 0.043043, 0.042872, 0.042932, 0.043291, 0.044029, 0.045101, 0.04658,
2023 0.048452, 0.050507, 0.052293, 0.053458, 0.054021, 0.053378, 0.052139, 0.53458, 0.050477, 0.048788, 0.047383,
2024 0.046341, 0.045631, 0.045178, 0.045022, 0.045112, 0.045395, 0.045833, 0.046402, 0.047084, 0.047865, 0.048726,
2025 0.049651, 0.050629, 0.051654, 0.052718, 0.053816, 0.054944, 0.056098, 0.057274, 0.058469, 0.059682, 0.060909,
2026 0.062149, 0.0634, 0.064661, 0.06593, 0.067207, 0.06849, 0.069778, 0.07107, 0.072367, 0.073666, 0.074968,
2027 0.076272, 0.077577, 0.078883, 0.080189, 0.081495, 0.082801, 0.084104, 0.085407, 0.086707, 0.088004, 0.089297,
2028 0.090585, 0.091867, 0.093142, 0.094408, 0.095664, 0.096907, 0.098134, 0.099336, 0.10051, 0.10164, 0.10273,
2029 0.10375, 0.10468, 0.10548, 0.10611, 0.10649, 0.10655, 0.10608, 0.10566, 0.1072, 0.10799, 0.10875,
2030 0.11103, 0.11491, 0.11819, 0.12051, 0.12211, 0.12339, 0.12449, 0.12556, 0.12663, 0.12771, 0.12881,
2031 0.12995, 0.13111, 0.13229, 0.13348, 0.13468, 0.13589, 0.13711, 0.13834, 0.13957, 0.1408, 0.14204,
2032 0.14328, 0.14452, 0.14576, 0.14701, 0.14825, 0.1495, 0.15075, 0.152, 0.15325, 0.1545, 0.15576,
2033 0.15701, 0.15826, 0.15952, 0.16077, 0.16203, 0.16328, 0.16454, 0.16579, 0.16705, 0.1683, 0.16956,
2034 0.17082, 0.17207, 0.17333, 0.17458, 0.17584, 0.1771, 0.17835, 0.17961, 0.18087, 0.18212, 0.18338,
2035 0.18464, 0.18589, 0.18715, 0.18841, 0.18966, 0.19092, 0.19218, 0.19343, 0.19469, 0.19595, 0.19721,
2036 0.19846, 0.19972, 0.20098, 0.20223, 0.20349, 0.20475, 0.20601, 0.20726, 0.20852, 0.20978, 0.21103,
2037 0.21229, 0.21355, 0.2148, 0.21606, 0.21732, 0.21857, 0.21983, 0.22109, 0.22234, 0.2236, 0.22486,
2038 0.22612, 0.22737, 0.22863, 0.22989, 0.23114, 0.2324, 0.23366, 0.23491, 0.23617, 0.2363};
2039
2040
2041 // Int_t xindex=0;
2042
2043 return drifttime[Int_t(x)];
2044}
2045
2046
2047//_______________________________________________________________________________
2048Double_t CbmTrdModuleSimR::GetStep(Double_t dist, Int_t roll)
2049{
2050 Double_t prob = 0.;
2051 Int_t steps = 1000;
2052 Double_t CalcGamma = 0.;
2053
2054 std::pair<Double_t, Double_t> bethe[12] = {
2055 std::make_pair(1.5, 1.5), std::make_pair(2, 1.1), std::make_pair(3, 1.025), std::make_pair(4, 1),
2056 std::make_pair(10, 1.1), std::make_pair(20, 1.2), std::make_pair(100, 1.5), std::make_pair(200, 1.6),
2057 std::make_pair(300, 1.65), std::make_pair(400, 1.675), std::make_pair(500, 1.7), std::make_pair(1000, 1.725)};
2058
2059 for (Int_t n = 0; n < 12; n++) {
2060 if (fGamma < bethe[0].first) {
2061 CalcGamma = bethe[0].second;
2062 break;
2063 }
2064 if (n == 11) {
2065 CalcGamma = bethe[11].second;
2066 break;
2067 }
2068
2069 if (fGamma >= bethe[n].first && fGamma <= bethe[n + 1].first) {
2070 Double_t dx = bethe[n + 1].first - bethe[n].first;
2071 Double_t dy = bethe[n + 1].second - bethe[n].second;
2072 Double_t slope = dy / dx;
2073 CalcGamma = (fGamma - bethe[n].first) * slope + bethe[n].second;
2074 break;
2075 }
2076 }
2077
2078 Double_t s = 0;
2079 Double_t D = 1 / (20.5 * CalcGamma);
2080 for (Int_t i = 1; i < steps; i++) {
2081 s = (dist / steps) * i;
2082 prob = (1 - TMath::Exp(-s / D)) * 100;
2083 if (prob >= roll) return s;
2084 }
2085
2086 return s;
2087}
2088
2089std::pair<Int_t, std::vector<Double_t>> CbmTrdModuleSimR::GetTotalSteps(Double_t In[3], Double_t Out[3], Double_t dist)
2090{
2091 Double_t prob = 0.;
2092 Int_t steps = 1000;
2093 Double_t CalcGamma = 0.;
2094 Double_t roll = gRandom->Integer(100);
2095
2096 std::pair<Double_t, Double_t> bethe[12] = {
2097 std::make_pair(1.5, 1.5), std::make_pair(2, 1.1), std::make_pair(3, 1.025), std::make_pair(4, 1),
2098 std::make_pair(10, 1.1), std::make_pair(20, 1.2), std::make_pair(100, 1.5), std::make_pair(200, 1.6),
2099 std::make_pair(300, 1.65), std::make_pair(400, 1.675), std::make_pair(500, 1.7), std::make_pair(1000, 1.725)};
2100
2101 for (Int_t n = 0; n < 12; n++) {
2102 if (fGamma < bethe[0].first) {
2103 CalcGamma = bethe[0].second;
2104 break;
2105 }
2106 if (n == 11) {
2107 CalcGamma = bethe[11].second;
2108 break;
2109 }
2110
2111 if (fGamma >= bethe[n].first && fGamma <= bethe[n + 1].first) {
2112 Double_t dx = bethe[n + 1].first - bethe[n].first;
2113 Double_t dy = bethe[n + 1].second - bethe[n].second;
2114 Double_t slope = dy / dx;
2115 CalcGamma = (fGamma - bethe[n].first) * slope + bethe[n].second;
2116 break;
2117 }
2118 }
2119
2120 Double_t pos[3] = {In[0], In[1], In[2]};
2121 Double_t lastpos[3] = {In[0], In[1], In[2]};
2122 Int_t pointcount = 0;
2123 std::vector<Double_t> posvec;
2124 Double_t D = 1 / (20.5 * CalcGamma);
2125 Double_t delta[3];
2126 Double_t s = 0;
2127 for (Int_t i = 0; i < 3; i++)
2128 delta[i] = (Out[i] - In[i]);
2129
2130 roll = gRandom->Integer(100);
2131 for (Int_t i = 1; i < steps; i++) {
2132 s = (dist / steps) * i;
2133 prob = (1 - TMath::Exp(-s / D)) * 100;
2134 if (prob >= roll) {
2135 Double_t move = 2 * (s / dist);
2136 for (Int_t n = 0; n < 3; n++)
2137 pos[n] = lastpos[n] + move * delta[n];
2138 for (Int_t n = 0; n < 3; n++)
2139 lastpos[n] = pos[n];
2140 // Double_t r = TMath::Sqrt((In[0]-pos[0])*(In[0]-pos[0])+(In[1]-pos[1])*(In[1]-pos[1]));
2141 if (TMath::Abs(pos[0]) < fDigiPar->GetSizeX() && TMath::Abs(pos[1]) < fDigiPar->GetSizeY() && (pos[2]) < Out[2]) {
2142 posvec.push_back(pos[0]);
2143 posvec.push_back(pos[1]);
2144 posvec.push_back(pos[2]);
2145 pointcount++;
2146 i = 1;
2147 roll = gRandom->Integer(100);
2148 continue;
2149 }
2150 }
2151 if (TMath::Abs(pos[0]) < fDigiPar->GetSizeX() && TMath::Abs(pos[1]) < fDigiPar->GetSizeY() && (pos[2]) < Out[2]) {
2152 continue;
2153 }
2154 break;
2155 }
2156
2157 return std::make_pair(pointcount, posvec);
2158}
2159
2160
TClonesArray * points
ClassImp(CbmConverterManager)
Helper class to convert unique channel ID back and forth.
TRD digitizer. Updated 24/04/2013 by Andrey Lebedev andrey.lebedev@gsi.de Updated 4/06/2018 by Alex B...
friend fvec sqrt(const fvec &a)
static constexpr size_t size()
Definition KfSimdPseudo.h:2
fscal v[fmask::Size]
Definition KfSimdPseudo.h:4
bool first
float Float_t
int Int_t
bool Bool_t
Generates beam ions for transport simulation.
int32_t GetNofLinks() const
Definition CbmMatch.h:42
void AddLink(const CbmLink &newLink)
Definition CbmMatch.cxx:47
Bookkeeping of time-slice content.
static uint32_t GetModuleId(uint32_t address)
Return module ID from address.
static uint32_t GetSectorId(uint32_t address)
Return sector ID from address.
static uint32_t GetColumnId(uint32_t address)
Return column ID from address.
static uint32_t GetLayerId(uint32_t address)
Return layer ID from address.
static uint32_t GetAddress(int32_t layerId, int32_t moduleId, int32_t sectorId, int32_t rowId, int32_t columnId)
Return address from system ID, layer, module, sector, column and row IDs.
static uint32_t GetRowId(uint32_t address)
Return row ID from address.
void SetTriggerType(const eTriggerType triggerType)
Set digi trigger type.
void SetFlag(const int32_t iflag, bool set=true)
Generic flag status setter.
static float Clk(eCbmTrdAsicType ty)
DAQ clock accessor for each ASIC.
Definition CbmTrdDigi.h:110
void SetErrorClass(const int32_t n)
Set digi error class (SPADIC only)
Definition CbmTrdDigi.h:235
double GetTime() const
Getter for physical time [ns]. Accounts for clock representation of each ASIC. In SPADIC case physica...
Definition CbmTrdDigi.h:155
double GetCharge() const
Common purpose charge getter.
void SetCharge(float c)
Charge setter for SPADIC ASIC.
static Bool_t AddNoise()
static Bool_t UseWeightedDist()
static Bool_t IsTimeBased()
Char_t fLayerId
layer identifier
const CbmTrdParModDigi * fDigiPar
read-out description of module
virtual Double_t GetDx() const
Shortcut getter size x/2 [cm].
virtual const Char_t * GetPath() const
UShort_t fModAddress
unique identifier for current module
CbmTrdParModAsic * fAsicPar
the set of ASIC operating on the module (owned)
virtual Double_t GetDy() const
Shortcut getter size y/2 [cm].
Simulation module implementation for rectangular pad geometry.
void SetPulsePars(Int_t mode)
std::pair< Int_t, std::vector< Double_t > > GetTotalSteps(Double_t In[3], Double_t Out[3], Double_t dist)
Int_t FlushBuffer(ULong64_t time=0)
Flush local digi buffer.
void SetSpadicResponse(Double_t calibration, Double_t tau)
void SetNoiseLevel(Double_t sigma_keV)
Bool_t MakeDigi(CbmTrdPoint *p, Double_t time, Bool_t TR)
Steering routine for converting MC point to digits.
void AddDigi(Int_t address, Double_t charge, Double_t chargeTR, Double_t time, Int_t trigger)
Double_t AddNoise(Double_t charge)
void CleanUp(Bool_t EB)
Double_t CalcPRF(Double_t x, Double_t W, Double_t h)
Double_t CalcResponse(Double_t t)
void AddDigitoBuffer(Int_t address, Double_t charge, Double_t chargeTR, Double_t time, Int_t trigger)
void ProcessBuffer(Int_t address)
std::map< Int_t, std::vector< std::vector< Int_t > > > fMCBuffer
void NoiseTime(ULong64_t eventTime)
CbmTrdCheckUtil * fQA
void ScanPadPlane(const Double_t *local_point, Double_t reldrift, Double_t clusterELoss, Double_t clusterELossTR, Int_t epoints, Int_t ipoint)
std::vector< Double_t > AddCorrelatedNoise(std::vector< Double_t > pulse)
CbmTrdRawToDigiR * fMessageConverter
CbmTrdModuleSimR(Int_t mod, Int_t ly, Int_t rot)
std::map< Int_t, std::vector< std::pair< CbmTrdDigi *, CbmMatch * > > > fAnalogBuffer
CbmTimeSlice * fTimeSlice
link to CBM time slice
Int_t CheckTrigger(std::vector< Double_t > pulse)
void AddDigitoPulseBuffer(Int_t address, Double_t reldrift, Double_t charge, Double_t chargeTR, Double_t time, Int_t trigger, Int_t epoints, Int_t ipoint)
std::map< Int_t, Double_t > fShiftQA
std::map< Int_t, Double_t > fMCQA
Int_t AddCrosstalk(Double_t address, Int_t i, Int_t sec, Int_t row, Int_t col, Int_t ncols)
void SetAsicPar(CbmTrdParModAsic *p=NULL)
std::vector< Double_t > MakePulse(Double_t charge, std::vector< Double_t > pulse, Int_t address)
Double_t AddDrifttime(Double_t x, Double_t z)
std::map< Int_t, std::vector< std::map< TString, Int_t > > > fLinkQA
void CheckTime(Int_t address)
Int_t GetMultiBin(std::vector< Double_t > pulse)
Bool_t CheckMulti(Int_t address, std::vector< Double_t > pulse)
std::map< Int_t, std::pair< Double_t, Int_t > > fMultiBuffer
std::map< Int_t, std::pair< std::vector< Double_t >, CbmMatch * > > fPulseBuffer
std::map< Int_t, Double_t > fTimeBuffer
void ProcessPulseBuffer(Int_t address, Bool_t FNcall, Bool_t MultiCall, Bool_t down, Bool_t up)
void SetPadPlaneScanArea(Int_t row)
void SetDistributionPoints(Int_t points)
void SetPulseMode(Bool_t pulsed)
Double_t GetStep(Double_t dist, Int_t roll)
void AddToPulse(Int_t address, Double_t charge, Double_t reldrift, std::vector< Double_t > pulse)
void CheckBuffer(Bool_t EB)
Bool_t DistributeCharge(Double_t pointin[3], Double_t pointout[3], Double_t delta[3], Double_t pos[3], Int_t ipoints)
std::shared_ptr< CbmTrdRadiator > fRadiator
Pointer to digitizer.
CbmTrdModuleSim()
Default constructor.
Int_t fInputId
MC input file number.
Int_t fPointId
MC point id being processed.
std::map< Int_t, std::pair< CbmTrdDigi *, CbmMatch * > > fDigiMap
Temporary storage for complete digis for each CBM address.
virtual void SetPositionMC(Double_t pos[3])
Int_t fEventId
MC event id being processed.
Double_t fXYZ[3]
MC position of the point in module coordinates.
virtual void SetChannelAddress(Int_t address)
Describe TRD module ASIC settings (electronic gain, delays, etc)
Definition of SPADIC parameters.
double GetYOut() const
Definition CbmTrdPoint.h:69
double GetXIn() const
Definition CbmTrdPoint.h:65
double GetZIn() const
Definition CbmTrdPoint.h:67
double GetXOut() const
Definition CbmTrdPoint.h:68
double GetYIn() const
Definition CbmTrdPoint.h:66
double GetZOut() const
Definition CbmTrdPoint.h:70