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