CbmRoot
CbmTofDigitize.cxx
Go to the documentation of this file.
1 
6 #include "CbmTofDigitize.h"
7 
8 // C++ includes
9 #include <algorithm>
10 #include <cassert>
11 #include <iomanip>
12 
13 // ROOT includes
14 #include "TClonesArray.h"
15 #include "TDirectory.h"
16 #include "TF2.h"
17 #include "TGeoManager.h"
18 #include "TLine.h"
19 #include "TMath.h"
20 #include "TROOT.h"
21 #include "TRandom3.h"
22 #include "TVector3.h"
23 #include <TFile.h>
24 
25 // FAIR classes and includes
26 #include "FairEventHeader.h"
27 #include "FairLogger.h"
28 #include "FairMCEventHeader.h"
29 #include "FairRootManager.h"
30 #include "FairRunAna.h"
31 #include "FairRunSim.h"
32 #include "FairRuntimeDb.h"
33 
34 // CBM includes
35 #include "CbmLink.h"
36 #include "CbmMCTrack.h"
37 #include "CbmMatch.h"
38 #include "CbmTofAddress.h" // in cbmdata/tof
39 #include "CbmTofDigi.h" // in cbmdata/tof
40 #include "CbmTofPoint.h" // in cbmdata/tof
41 
42 // TOF includes
43 #include "CbmTofCell.h" // in tof/TofData
44 #include "CbmTofCreateDigiPar.h" // in tof/TofTools
45 #include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
46 #include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
47 #include "CbmTofDigiBdfPar.h" // in tof/TofParam
48 #include "CbmTofDigiPar.h" // in tof/TofParam
49 #include "CbmTofGeoHandler.h" // in tof/TofTools
50 
51 using std::cout;
52 using std::endl;
53 using std::fixed;
54 using std::right;
55 using std::setprecision;
56 using std::setw;
57 
58 // Gauss Integration Constants
59 const Int_t kiNbIntPts = 2;
60 Bool_t bFakeBeamCounter = kFALSE;
61 
62 /************************************************************************************/
63 struct CompTimesExp {
64  typedef std::pair<CbmTofDigi*, CbmMatch*> data;
65  bool operator()(data pair1, data pair2) {
66  Bool_t bOK = (pair1.first->GetTime() < pair2.first->GetTime());
67  /*
68  if(!bOK) {
69  LOG(debug) << Form(" ReSort digis 0x%08x , %f, %f", pDigi1->GetAddress(), pDigi1->GetTime(), pDigi2->GetTime() )
70  ;
71  }
72  */
73  return bOK;
74  }
75 } CompTExp;
76 
77 // If enabled add an offset depending on the strip length by making the full propagation time
78 // to each strips end
79 // If disabled just the time to the center is used
80 // TODO: Add this as parameter
81 //#define FULL_PROPAGATION_TIME
82 Int_t fCurrentMCEntry;
83 
84 /************************************************************************************/
85 CbmTofDigitize::CbmTofDigitize()
86  : CbmDigitize<CbmTofDigi>("TofDigitize")
87  , fdFeeGainSigma(0.)
88  , fdFeeTotThr(0.)
89  , fdTimeResElec(0.)
90  , fdSignalPropSpeed(0.)
91  , fh1ClusterSizeProb()
92  , fh1ClusterTotProb()
93  , fvdSignalVelocityRpc()
94  , fdChannelGain()
95  , fGeoHandler(new CbmTofGeoHandler())
96  , fTofId(NULL)
97  , fDigiPar(NULL)
98  , fChannelInfo(NULL)
99  , fDigiBdfPar(NULL)
100  , fiNbElecChTot(0)
101  , fvRpcChOffs()
102  , fTofPointsColl(NULL)
103  , fMcTracksColl(NULL)
104  , fStorDigi()
105  , fStorDigiMatch()
106  , fvlTrckChAddr()
107  , fvlTrckRpcAddr()
108  , fvlTrckRpcTime()
109  , fiNbDigis(0)
110  , fsHistoOutFilename("")
111  , fhTofPointsPerTrack(NULL)
112  , fhTofPtsInTrkVsGapInd(NULL)
113  , fhTofPtsInTrkVsGapIndPrm(NULL)
114  , fhTofPtsInTrkVsGapIndSec(NULL)
115  , fhTofPtsPosVsGap()
116  , fhEvtProcTime(NULL)
117  , fhDigiMergeTime(NULL)
118  , fhDigiNbElecCh(NULL)
119  , fhProcTimeEvtSize(NULL)
120  , fhMeanDigiPerTrack(NULL)
121  , fhMeanFiredPerTrack(NULL)
122  , fhPtTime(NULL)
123  , fhDigiTime(NULL)
124  , fhDigiTimeRes(NULL)
125  , fhDigiTimeResB(NULL)
126  , fhToTDist(NULL)
127  , fhElecChOccup(NULL)
128  , fhMultiDigiEvtElCh(NULL)
129  , fhNbDigiEvtElCh(NULL)
130  , fhNbTracksEvtElCh(NULL)
131  , fhFiredEvtElCh(NULL)
132  , fhMultiProbElCh(NULL)
133  , fStart()
134  , fStop()
135  , fdDigitizeTime(0.)
136  , fdMergeTime(0.)
137  , fTimer()
138  , fiNofEvents(0)
139  , fdNofTofMcTrkTot(0.)
140  , fdNofPointsTot(0.)
141  , fdNofDigisTot(0.)
142  , fdTimeTot(0.)
143  , fsBeamInputFile("")
144  , fbMonitorHistos(kTRUE)
145  , fbMcTrkMonitor(kFALSE)
146  , fbTimeBasedOutput(kFALSE)
147  , fbAllowPointsWithoutTrack(kFALSE)
148  ,
149  //fiCurrentFileId(-1),
150  //fiCurrentEventId(-1),
151  //fdCurrentEventTime(0.),
152  fdDigiTimeConvFactor(1.) {}
153 
154 CbmTofDigitize::CbmTofDigitize(const char* name, Int_t)
155  : CbmDigitize<CbmTofDigi>(name)
156  , fdFeeGainSigma(0.)
157  , fdFeeTotThr(0.)
158  , fdTimeResElec(0.)
159  , fdSignalPropSpeed(0.)
160  , fh1ClusterSizeProb()
161  , fh1ClusterTotProb()
162  , fvdSignalVelocityRpc()
163  , fdChannelGain()
164  , fGeoHandler(new CbmTofGeoHandler())
165  , fTofId(NULL)
166  , fDigiPar(NULL)
167  , fChannelInfo(NULL)
168  , fDigiBdfPar(NULL)
169  , fiNbElecChTot(0)
170  , fvRpcChOffs()
171  , fTofPointsColl(NULL)
172  , fMcTracksColl(NULL)
173  , fStorDigi()
174  , fStorDigiMatch()
175  , fvlTrckChAddr()
176  , fvlTrckRpcAddr()
177  , fvlTrckRpcTime()
178  , fiNbDigis(0)
179  , fsHistoOutFilename("")
180  , fhTofPointsPerTrack(NULL)
181  , fhTofPtsInTrkVsGapInd(NULL)
182  , fhTofPtsInTrkVsGapIndPrm(NULL)
183  , fhTofPtsInTrkVsGapIndSec(NULL)
184  , fhTofPtsPosVsGap()
185  , fhEvtProcTime(NULL)
186  , fhDigiMergeTime(NULL)
187  , fhDigiNbElecCh(NULL)
188  , fhProcTimeEvtSize(NULL)
189  , fhMeanDigiPerTrack(NULL)
190  , fhMeanFiredPerTrack(NULL)
191  , fhPtTime(NULL)
192  , fhDigiTime(NULL)
193  , fhDigiTimeRes(NULL)
194  , fhDigiTimeResB(NULL)
195  , fhToTDist(NULL)
196  , fhElecChOccup(NULL)
197  , fhMultiDigiEvtElCh(NULL)
198  , fhNbDigiEvtElCh(NULL)
199  , fhNbTracksEvtElCh(NULL)
200  , fhFiredEvtElCh(NULL)
201  , fhMultiProbElCh(NULL)
202  , fStart()
203  , fStop()
204  , fdDigitizeTime(0.)
205  , fdMergeTime(0.)
206  , fTimer()
207  , fiNofEvents(0)
208  , fdNofTofMcTrkTot(0.)
209  , fdNofPointsTot(0.)
210  , fdNofDigisTot(0.)
211  , fdTimeTot(0.)
212  , fsBeamInputFile("")
213  , fbMonitorHistos(kTRUE)
214  , fbMcTrkMonitor(kFALSE)
215  , fbTimeBasedOutput(kFALSE)
216  , fbAllowPointsWithoutTrack(kFALSE)
217  ,
218  //fiCurrentFileId(-1),
219  //fiCurrentEventId(-1),
220  //fdCurrentEventTime(0.),
221  fdDigiTimeConvFactor(1.) {}
222 
223 CbmTofDigitize::~CbmTofDigitize() {
224  if (fGeoHandler) delete fGeoHandler;
225  // DeleteHistos(); // <-- if needed ?
226 }
227 
228 /************************************************************************************/
229 // FairTasks inherited functions
230 InitStatus CbmTofDigitize::Init() {
231  std::cout << std::endl;
232  LOG(info) << "==========================================================";
233  LOG(info) << GetName() << ": Initialisation";
234  if (fEventMode) LOG(info) << GetName() << ": Using event mode.";
235 
236  // If input file was not set explicitly, use default one
237  if (fsBeamInputFile.IsNull()) {
238  TString fileName = gSystem->Getenv("VMCWORKDIR");
239  fileName += "/parameters/tof/test_bdf_input.root";
240  SetInputFileName(fileName);
241  LOG(info) << GetName() << ": Using default parameter file " << fileName;
242  }
243 
244  if (kFALSE == RegisterInputs()) return kFATAL;
245 
246  RegisterOutput();
247 
248  if (kFALSE == InitParameters()) return kFATAL;
249 
250  if (kFALSE == LoadBeamtimeValues()) return kFATAL;
251 
252  if (kFALSE == CreateHistos()) return kFATAL;
253 
254  LOG(info) << GetName() << ": Initialisation successful";
255  LOG(info) << "==========================================================";
256  return kSUCCESS;
257 }
258 
259 void CbmTofDigitize::SetParContainers() {
260  LOG(info) << " CbmTofDigitize => Get the digi parameters for tof";
261 
262  // Get Base Container
263  FairRun* ana = FairRun::Instance();
264  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
265 
266  fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
267 
268  fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
269 }
270 
271 void CbmTofDigitize::Exec(Option_t* /*option*/) {
272  // --- Start timer and reset counters
273  fTimer.Start();
274 
275  // Get FairEventHeader information for CbmLink objects
276  GetEventInfo();
277  LOG(debug) << fName << ": Processing event " << fCurrentEvent << " at "
278  << fCurrentEventTime << " ns";
279 
280 
281  fiNbDigis = 0;
282  Int_t nPoints = fTofPointsColl->GetEntriesFast();
283 
284  fStart.Set();
285  LOG(debug) << fName << ": Using cluster model "
286  << fDigiBdfPar->GetClusterModel();
287  switch (fDigiBdfPar->GetClusterModel()) {
288  case 0: DigitizeDirectClusterSize(); break;
289  case 1: DigitizeFlatDisc(); break;
290  case 2: DigitizeGaussCharge(); break;
291  default: DigitizeDirectClusterSize(); break;
292  }
293 
294  fStop.Set();
295  fdDigitizeTime = fStop.GetSec() - fStart.GetSec()
296  + (fStop.GetNanoSec() - fStart.GetNanoSec()) / 1e9;
297 
298  LOG(debug) << fName << ": Merge digis ";
299  MergeSameChanDigis();
300  fStart.Set();
301  fdMergeTime = fStart.GetSec() - fStop.GetSec()
302  + (fStart.GetNanoSec() - fStop.GetNanoSec()) / 1e9;
303 
304  // --- Counters
305  fTimer.Stop();
306  fiNofEvents++;
307  fdNofPointsTot += nPoints;
308  fdNofDigisTot += fiNbDigis;
309  fdTimeTot += fTimer.RealTime();
310 
311  // --- Event log
312  LOG(info) << "+ " << setw(15) << GetName() << ": Event " << setw(6) << right
313  << fCurrentEvent << " at " << fixed << setprecision(3)
314  << fCurrentEventTime << " ns, points: " << nPoints
315  << ", digis: " << fiNbDigis << ". Exec time " << setprecision(6)
316  << fTimer.RealTime() << " s.";
317 }
318 
319 void CbmTofDigitize::Finish() {
320  std::cout << std::endl;
321  LOG(info) << "=====================================";
322 
323  // Final printout for reference
324  LOG(info) << GetName() << ": Run summary (Time includes Hist filling)";
325  LOG(info) << "Events processed : " << fiNofEvents;
326  if (fbMcTrkMonitor) {
327  LOG(info) << "Tof MC Track / event: " << fdNofPointsTot / fdNofTofMcTrkTot;
328  } // if( fbMcTrkMonitor )
329  LOG(info) << "TofPoint / event : "
330  << fdNofPointsTot / static_cast<Double_t>(fiNofEvents);
331  LOG(info) << "TofDigi / event : "
332  << fdNofDigisTot / static_cast<Double_t>(fiNofEvents);
333  LOG(info) << "Digis per point : " << fdNofDigisTot / fdNofPointsTot;
334  if (fbMcTrkMonitor) {
335  LOG(info) << "Digis per TofMcTrack: " << fdNofDigisTot / fdNofTofMcTrkTot;
336  } // if( fbMcTrkMonitor )
337  LOG(info) << "Real time per event : "
338  << fdTimeTot / static_cast<Double_t>(fiNofEvents) << " s";
339 
340  WriteHistos();
341  // Prevent them from being sucked in by the CbmHadronAnalysis WriteHistograms method
342  DeleteHistos();
343 
344  // Clear the cluster size and efficiency histos copies inside the parameter
345  fDigiBdfPar->ClearHistos();
346  LOG(info) << "=====================================";
347 }
348 
349 /************************************************************************************/
350 // Functions common for all clusters approximations
351 Bool_t CbmTofDigitize::RegisterInputs() {
352  FairRootManager* fManager = FairRootManager::Instance();
353 
354  fTofPointsColl = (TClonesArray*) fManager->GetObject("TofPoint");
355  if (NULL == fTofPointsColl) {
356  LOG(error) << "CbmTofDigitize::RegisterInputs => Could not get the "
357  "TofPoint TClonesArray!!!";
358  return kFALSE;
359  } // if( NULL == fTofPointsColl)
360 
361  fMcTracksColl = (TClonesArray*) fManager->GetObject("MCTrack");
362  if (NULL == fMcTracksColl) {
363  LOG(error) << "CbmTofDigitize::RegisterInputs => Could not get the MCTrack "
364  "TClonesArray!!!";
365  return kFALSE;
366  } // if( NULL == fMcTracksColl)
367 
368  return kTRUE;
369 }
370 Bool_t CbmTofDigitize::InitParameters() {
371  // Initialize the TOF GeoHandler
372  Bool_t isSimulation = kFALSE;
373  Int_t iGeoVersion = fGeoHandler->Init(isSimulation);
374  if (k12b > iGeoVersion) {
375  LOG(error) << "CbmTofDigitize::InitParameters => Only compatible with "
376  "geometries after v12b !!!";
377  return kFALSE;
378  }
379 
380  LOG(info) << fName << ": GeoVersion " << iGeoVersion;
381 
382  switch (iGeoVersion) {
383  case k12b: fTofId = new CbmTofDetectorId_v12b(); break;
384  case k14a: fTofId = new CbmTofDetectorId_v14a(); break;
385  default:
386  LOG(error) << "CbmTofDigitize::InitParameters: Invalid Detector ID "
387  << iGeoVersion;
388  }
389 
390  // create digitization parameters from geometry file
391  CbmTofCreateDigiPar* tofDigiPar =
392  new CbmTofCreateDigiPar("TOF Digi Producer", "TOF task");
393  LOG(info) << "Create DigiPar";
394  tofDigiPar->Init();
395 
396  // Printout option for crosscheck
397  LOG(debug) << " CbmTofDigitize::LoadBeamtimeValues Digi Par contains "
398  << fDigiPar->GetNrOfModules() << " cells ";
399  return kTRUE;
400 }
401 Bool_t CbmTofDigitize::LoadBeamtimeValues() {
402  LOG(info) << fName << ": Load beamtime values from " << fsBeamInputFile;
403 
404  fDigiBdfPar->SetInputFile(fsBeamInputFile);
405 
406  // Add Param printout only if DEBUG level ON
407  if (gLogger->IsLogNeeded(fair::Severity::debug)) fDigiBdfPar->printParams();
408 
409  if (kFALSE == fDigiBdfPar->LoadBeamtimeHistos()) {
410  LOG(fatal) << "CbmTofDigitize::LoadBeamtimeValues: Cluster properties from "
411  "testbeam could not be loaded! ";
412  return kFALSE;
413  }
414 
415  // // Printout option for crosscheck
416  // fDigiBdfPar->printParams();
417 
418  // Obtain some of the constants
419  fdFeeGainSigma = fDigiBdfPar->GetFeeGainSigma();
420  fdFeeTotThr = fDigiBdfPar->GetFeeThreshold();
421  fdTimeResElec = fDigiBdfPar->GetFeeTimeRes();
422  fdSignalPropSpeed = fDigiBdfPar->GetSignalSpeed();
423 
424  // Generate the gain/fee channel matrix
425  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
426 
427  LOG(debug2)
428  << "CbmTofDigitize::LoadBeamtimeValues: ini gain values for SmTypes "
429  << iNbSmTypes;
430 
431  fdChannelGain.resize(iNbSmTypes);
432 
433  TRandom3 randFeeGain;
434  // randFeeGain.SetSeed(0);
435 
436  // Save the total number of electronic channels for the "Occupancy" estimation
437  // and the first channel of each RPC for the multiple signals estimation
438  fiNbElecChTot = 0;
439  fvRpcChOffs.resize(iNbSmTypes);
440  fvdSignalVelocityRpc.resize(iNbSmTypes);
441 
442  for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
443  LOG(debug2) << "CbmTofDigitize::LoadBeamtimeValues => Gain for SM type "
444  << iSmType;
445  Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
446  Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
447 
448  fvdSignalVelocityRpc[iSmType].resize(iNbSm);
449  fdChannelGain[iSmType].resize(iNbSm * iNbRpc);
450  fvRpcChOffs[iSmType].resize(iNbSm);
451 
452  if (iSmType == 5 && iNbSm > 0) {
453  bFakeBeamCounter = kTRUE;
454  LOG(info) << "Generate Fake Beam Counter digis";
455  }
456 
457  for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
458  fvdSignalVelocityRpc[iSmType][iSm].resize(iNbRpc);
459  for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++)
460  if (0.0 < fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc))
461  fvdSignalVelocityRpc[iSmType][iSm][iRpc] = fDigiBdfPar->GetSigVel(
462  iSmType, iSm, iRpc); // convert into cm/ns if necessary (FIXME)
463  else
464  fvdSignalVelocityRpc[iSmType][iSm][iRpc] = fdSignalPropSpeed;
465  }
466 
467  for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
468  fvRpcChOffs[iSmType][iSm].resize(iNbRpc);
469 
470  for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
471  LOG(debug2) << "CbmTofDigitize::LoadBeamtimeValues => Gain for SM/Rpc "
472  << iSm << "/" << iRpc;
473  Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
474  Int_t iNbSides = 2 - fDigiBdfPar->GetChanType(iSmType, iRpc);
475 
476  // Save the total number of electronic channels for the "Occupancy" estimation
477  // and the first channel of each RPC for the multiple signals estimation
478  fiNbElecChTot += iNbCh * iNbSides;
479  fvRpcChOffs[iSmType][iSm][iRpc] = fiNbElecChTot;
480 
481  fdChannelGain[iSmType][iSm * iNbRpc + iRpc].resize(iNbCh * iNbSides);
482 
483  for (Int_t iCh = 0; iCh < iNbCh; iCh++)
484  for (Int_t iSide = 0; iSide < iNbSides; iSide++) {
485  if (0 < fdFeeGainSigma) {
486  fdChannelGain[iSmType][iSm * iNbRpc + iRpc]
487  [iCh * iNbSides + iSide] =
488  randFeeGain.Gaus(1.0, fdFeeGainSigma);
489  if (fdChannelGain[iSmType][iSm * iNbRpc + iRpc]
490  [iCh * iNbSides + iSide]
491  < 0.0)
492  LOG(error) << "CbmTofDigitize::LoadBeamtimeValues => Neg. Gain "
493  "for SM/Rpc "
494  << iSm << "/" << iRpc << " "
495  << fdChannelGain[iSmType][iSm * iNbRpc + iRpc]
496  [iCh * iNbSides + iSide];
497  } // if( 0 < fdFeeGainSigma )
498  else
499  fdChannelGain[iSmType][iSm * iNbRpc + iRpc]
500  [iCh * iNbSides + iSide] = 1;
501  //if(iSmType==5) fdChannelGain[iSmType][iSm*iNbRpc + iRpc][iCh*iNbSides + iSide] *= 10.; //FIXME, Diamond
502  }
503  } // for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ )
504  } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ )
505  } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
506 
507  TDirectory* oldir = gDirectory;
508  gROOT->cd();
509  // Generate the Probability conversion histograms for Cluster size
510  // and cluster charge
511  fh1ClusterSizeProb.resize(iNbSmTypes);
512  fh1ClusterTotProb.resize(iNbSmTypes);
513  for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
514  LOG(debug2) << "CbmTofDigitize::LoadBeamtimeValues => SM type " << iSmType;
515 
516  if (0 == fDigiBdfPar->GetClusterModel()) {
517  // Generating the Cluster Radius Probability from the Cluster Size distribution
518  fh1ClusterSizeProb[iSmType] =
519  new TH1D(Form("ClusterSizeProb%03d", iSmType),
520  "Random throw to Cluster size mapping",
521  10000,
522  0.0,
523  1.0);
524  TH1* hClusterSize = fDigiBdfPar->GetClustSizeHist(iSmType);
525  Int_t iNbBinsSize = hClusterSize->GetNbinsX();
526  Double_t dNbBinsProb = fh1ClusterSizeProb[iSmType]->GetNbinsX();
527  Int_t iProbBin = 1;
528  Double_t dIntegral = 0.;
529  // Double_t dIntSize = hClusterSize->GetEntries();
530  // BUGFIX
531  // The histogram integral is the sum of bin contents (1-N), while the
532  // number of histogram entries is given by the number of times the
533  // Fill() method was called.
534  // The quantile histogram 'fh1ClusterSizeProb' computed below might
535  // contain a few F^-1(u) = 0 towards u = 1 if the underflow and/or
536  // overflow bins of the original PDF histogram were populated. This
537  // in return would lead to an excess at 0 in the spectrum of values
538  // sampled from the quantile histogram.
539  Double_t dIntSize = hClusterSize->Integral();
540  Double_t dSizeVal = 0.;
541 
542  for (Int_t iBin = 1; iBin <= iNbBinsSize; iBin++) {
543  dIntegral += hClusterSize->GetBinContent(iBin);
544  if ((Double_t) iProbBin / dNbBinsProb < dIntegral / dIntSize)
545  dSizeVal = hClusterSize->GetBinCenter(iBin);
546  while ((Double_t) iProbBin / dNbBinsProb < dIntegral / dIntSize) {
547  fh1ClusterSizeProb[iSmType]->SetBinContent(iProbBin, dSizeVal);
548  iProbBin++;
549  } // while( (Double_t)iProbBin/dNbBinsProb < dIntegral/dIntTot )
550  } // for( Int_t iBin = 1; iBin <= iNbBinsTot; iBin++ )
551  fh1ClusterSizeProb[iSmType]->SetBinContent(iProbBin, dSizeVal);
552  } // if( 0 == fDigiBdfPar->GetClusterModel() )
553  else
554  switch (fDigiBdfPar->GetClusterRadiusModel()) {
555  case 0: {
556  // Single value using the beamtime cluster size distribution mean
557  // => Just copy pointer to the beamtime histogram
558  fh1ClusterSizeProb[iSmType] = fDigiBdfPar->GetClustSizeHist(iSmType);
559  LOG(debug) << "CbmTofDigitize::LoadBeamtimeValues => SM type "
560  << iSmType << " Mean Cluster Size "
561  << (fh1ClusterSizeProb[iSmType]->GetMean());
562  break;
563  }
564  case 1: {
565  // from beamtime cluster size distribution
566  // Ingo tip: use Landau distribution until it matchs
567  fh1ClusterSizeProb[iSmType] = fDigiBdfPar->GetClustSizeHist(iSmType);
568  LOG(debug) << "CbmTofDigitize::LoadBeamtimeValues => SM type "
569  << iSmType << " Mean Cluster Size "
570  << (fh1ClusterSizeProb[iSmType]->GetMean());
571  break;
572  }
573  default: {
574  // Should not happen !
575  fh1ClusterSizeProb[iSmType] = 0;
576  break;
577  }
578  } // switch( fDigiBdfPar->GetClusterRadiusModel() )
579 
580 
581  fh1ClusterTotProb[iSmType] = new TH1D(Form("ClusterTotProb%03d", iSmType),
582  "Random throw to Cluster Tot mapping",
583  10000,
584  0.0,
585  1.0);
586  TH1* hClusterTot = fDigiBdfPar->GetClustTotHist(iSmType);
587  Int_t iNbBinsTot = hClusterTot->GetNbinsX();
588  Double_t dNbBinsProb = fh1ClusterTotProb[iSmType]->GetNbinsX();
589  Int_t iProbBin = 1;
590  Double_t dIntegral = 0.;
591  // Double_t dIntTot = hClusterTot->GetEntries();
592  // BUGFIX
593  // The histogram integral is the sum of bin contents (1-N), while the
594  // number of histogram entries is given by the number of times the
595  // Fill() method was called.
596  // The quantile histogram 'fh1ClusterTotProb' computed below might
597  // contain a few F^-1(u) = 0 towards u = 1 if the underflow and/or
598  // overflow bins of the original PDF histogram were populated. This
599  // in return would lead to an excess at 0 in the spectrum of values
600  // sampled from the quantile histogram.
601  Double_t dIntTot = hClusterTot->Integral();
602  Double_t dTotVal = 0.;
603  Double_t dPsToNs =
604  1e3; // FIXME? Most histograms from analysis are in ps, simulation used ns !!
605 
606  for (Int_t iBin = 1; iBin <= iNbBinsTot; iBin++) {
607  dIntegral += hClusterTot->GetBinContent(iBin);
608  if ((Double_t) iProbBin / dNbBinsProb < dIntegral / dIntTot)
609  dTotVal = hClusterTot->GetBinLowEdge(iBin) / dPsToNs;
610  while ((Double_t) iProbBin / dNbBinsProb < dIntegral / dIntTot) {
611  fh1ClusterTotProb[iSmType]->SetBinContent(iProbBin, dTotVal);
612  iProbBin++;
613  } // while( (Double_t)iProbBin/dNbBinsProb < dIntegral/dIntTot )
614  } // for( Int_t iBin = 1; iBin <= iNbBinsTot; iBin++ )
615  fh1ClusterTotProb[iSmType]->SetBinContent(iProbBin, dTotVal);
616  } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
617  gDirectory->cd(oldir->GetPath());
618 
619  // Prepare the intermediate Digi storage
620  fStorDigi.resize(iNbSmTypes);
621  fStorDigiMatch.resize(iNbSmTypes);
622  for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
623  Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
624  Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
625  fStorDigi[iSmType].resize(iNbSm * iNbRpc);
626  fStorDigiMatch[iSmType].resize(iNbSm * iNbRpc);
627  for (Int_t iSm = 0; iSm < iNbSm; iSm++)
628  for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
629  Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
630  Int_t iNbSides = 2 - fDigiBdfPar->GetChanType(iSmType, iRpc);
631  fStorDigi[iSmType][iSm * iNbRpc + iRpc].resize(iNbCh * iNbSides);
632  fStorDigiMatch[iSmType][iSm * iNbRpc + iRpc].resize(iNbCh * iNbSides);
633  } // for each (Sm, rpc) pair
634  } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
635 
636  LOG(debug) << "CbmTofDigitize::LoadBeamtimeValues: GeoVersion "
637  << fGeoHandler->GetGeoVersion();
638 
639  // TEST stupid stuffs
640  Int_t iNbGeantChannels = fDigiPar->GetNrOfModules();
641  Int_t iMinSmType = 17; // 4b for sm type in ID v12b normally
642  Int_t iMaxSmType = -1; // 4b for sm type in ID v12b normally
643  Int_t iMinSm = 257; // 8b for sm type in ID v12b normally
644  Int_t iMaxSm = -1; // 8b for sm type in ID v12b normally
645  Int_t iMinRpc = 9; // 3b for sm type in ID v12b normally
646  Int_t iMaxRpc = -1; // 3b for sm type in ID v12b normally
647  Int_t iMinGap = 17; // 4b for sm type in ID v12b normally
648  Int_t iMaxGap = -1; // 4b for sm type in ID v12b normally
649  Int_t iMinCell = 257; // 8b for sm type in ID v12b normally
650  Int_t iMaxCell = -1; // 4b for sm type in ID v12b normally
651  for (Int_t iCellInd = 0; iCellInd < iNbGeantChannels; iCellInd++) {
652  Int_t iCellId = fDigiPar->GetCellId(iCellInd);
653 
654  if (iMinSmType > fTofId->GetSMType(iCellId))
655  iMinSmType = fTofId->GetSMType(iCellId);
656  if (iMaxSmType < fTofId->GetSMType(iCellId))
657  iMaxSmType = fTofId->GetSMType(iCellId);
658 
659  if (iMinSm > fTofId->GetSModule(iCellId))
660  iMinSm = fTofId->GetSModule(iCellId);
661  if (iMaxSm < fTofId->GetSModule(iCellId))
662  iMaxSm = fTofId->GetSModule(iCellId);
663 
664  if (iMinRpc > fTofId->GetCounter(iCellId))
665  iMinRpc = fTofId->GetCounter(iCellId);
666  if (iMaxRpc < fTofId->GetCounter(iCellId))
667  iMaxRpc = fTofId->GetCounter(iCellId);
668 
669  if (iMinGap > fTofId->GetGap(iCellId)) iMinGap = fTofId->GetGap(iCellId);
670  if (iMaxGap < fTofId->GetGap(iCellId)) iMaxGap = fTofId->GetGap(iCellId);
671 
672  if (iMinCell > fTofId->GetCell(iCellId))
673  iMinCell = fTofId->GetCell(iCellId);
674  if (iMaxCell < fTofId->GetCell(iCellId))
675  iMaxCell = fTofId->GetCell(iCellId);
676  }
677  LOG(debug) << "SM type min " << iMinSmType << " max " << iMaxSmType;
678  LOG(debug) << "SM min " << iMinSm << " max " << iMaxSm;
679  LOG(debug) << "Rpc min " << iMinRpc << " max " << iMaxRpc;
680  LOG(debug) << "Gap min " << iMinGap << " max " << iMaxGap;
681  LOG(debug) << "Chan min " << iMinCell << " max " << iMaxCell;
682  //
683  return kTRUE;
684 }
685 /************************************************************************************/
686 // Histogramming functions
687 Bool_t CbmTofDigitize::CreateHistos() {
688  if (!fbMonitorHistos) { return kTRUE; }
689 
690  TDirectory* oldir =
691  gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager!
692  gROOT
693  ->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !
694 
695  fhTofPointsPerTrack = new TH1I("TofDigiBdf_TofPtPerTrk",
696  "Number of Tof Points per MC track reaching "
697  "the TOF; Nb Tof Points in Tracks []",
698  16,
699  0.0,
700  16.0);
701  fhTofPtsInTrkVsGapInd =
702  new TH2I("TofDigiBdf_TofPtInTrkGap",
703  "Gap index vs Number of Tof Points in corresponding MC Track; Nb "
704  "Tof Pts in Track []; Gap Index []",
705  16,
706  0.0,
707  16.0,
708  10,
709  -0.5,
710  9.5);
711  fhTofPtsInTrkVsGapIndPrm =
712  new TH2I("TofDigiBdf_TofPtInTrkGapPrm",
713  "Gap index vs Number of Tof Points in corresponding MC Track, "
714  "primaries; Nb Tof Pts in Track []; Gap Index []",
715  16,
716  0.0,
717  16.0,
718  10,
719  -0.5,
720  9.5);
721  fhTofPtsInTrkVsGapIndSec =
722  new TH2I("TofDigiBdf_TofPtInTrkGapSec",
723  "Gap index vs Number of Tof Points in corresponding MC Track, "
724  "secondaries; Nb Tof Pts in Track []; Gap Index []",
725  16,
726  0.0,
727  16.0,
728  10,
729  -0.5,
730  9.5);
731  for (Int_t iGap = 0; iGap < 10; iGap++)
732  fhTofPtsPosVsGap[iGap] =
733  new TH2I(Form("TofDigiBdf_fhTofPtsPosVsGap%d", iGap),
734  Form("Tof Point positions in gap %d; X [cm]; Y [cm]", iGap),
735  750,
736  -750,
737  750,
738  500,
739  -500,
740  500);
741  /*
742  fhTofPointsPerTrackVsPdg = new TH2I( "TofDigiBdf_TofPtPerTrkVsPdg", "Number of Tof Points per MC track in each event; PDG []; Nb Points/Nb Tracks []",
743  800, -4000, 4000,
744  16, 0.0, 16.0 );
745  fhMeanPointPerTrack = new TH1I( "TofDigiBdf_PtPerTrk", "Mean Number of Tof Points per MC track in each event; Nb Points/Nb Tracks []",
746  800, 0.0, 16.0 );
747  fhMeanPointPerTrack2d = new TH2I( "TofDigiBdf_PtPerTrk2d", "Mean Number of Tof Points per MC track in each event; Nb Tracks []; Nb Tof Points/Nb Tracks []",
748  300, 0.0, 3000.0,
749  800, 0.0, 16.0 );
750  fhMeanDigiPerPoint = new TH1I( "TofDigiBdf_DigiPerPt", "Mean Number of ToF Digi per Tof Point in each event; Nb Digi/Nb Points []",
751  500, 0.0, 10.0 );
752  fhMeanFiredPerPoint = new TH1I( "TofDigiBdf_FirPerPt", "Mean Number of fired channel per Tof Point in each event; Nb Fired/Nb Points []",
753  500, 0.0, 10.0 );
754 */
755  fhEvtProcTime = new TH1I("TofDigiBdf_EvtProcTime",
756  "Time needed to process each event; Time [s]",
757  40000,
758  0.0,
759  40.0);
760  fhDigiMergeTime =
761  new TH1I("TofDigiBdf_DigiMergeTime",
762  "Time needed to merge the digis for each event; Time [s]",
763  4000,
764  0.0,
765  0.4);
766  fhDigiNbElecCh = new TH1I(
767  "TofDigiBdf_DigiNbElecCh",
768  "Nb of digis per channel before merging; Nb Digis in same elec. ch []",
769  50,
770  0.0,
771  50);
772  fhProcTimeEvtSize = new TH2I("TofDigiBdf_ProcTimeEvtSize",
773  "Time needed to process each event vs Event "
774  "Size; Event Size [TofPoints]; Time [s]",
775  600,
776  0.0,
777  12000.0,
778  400,
779  0.0,
780  4.0);
781  fhMeanDigiPerTrack = new TH1I(
782  "TofDigiBdf_DigiPerTrk",
783  "Mean Number of ToF Digi per MC track in each event; Nb Digi/Nb Tracks []",
784  500,
785  0.0,
786  10.0);
787  fhMeanFiredPerTrack = new TH1I("TofDigiBdf_FirPerTrk",
788  "Mean Number of fired channel per MC track in "
789  "each event; Nb Fired/Nb Tracks []",
790  500,
791  0.0,
792  10.0);
793  fhPtTime = new TH1I(
794  "TofDigiBdf_PtTime",
795  "TofPoints Time distribution, summed for all Pts/channels; Time [ns]",
796  10000,
797  0.0,
798  10000.0);
799  fhDigiTime =
800  new TH1I("TofDigiBdf_DigiTime",
801  "Time distribution, summed for all digis/channels; Time [ns]",
802  10000,
803  0.0,
804  10000.0);
805  fhDigiTimeRes = new TH1I("TofDigiBdf_DigiTimeRes",
806  "Time to True time distribution, summed for all "
807  "digis/channels; Time Digi - Time Pt [ns]",
808  10000,
809  -50.0,
810  50.0);
811  fhDigiTimeResB = new TH2I("TofDigiBdf_DigiTimeResB",
812  "Time to True time distribution, summed for all "
813  "digis/channels; Time Digi - Time Pt [ns]",
814  10000,
815  -50.0,
816  50.0,
817  6,
818  0,
819  6);
820  fhToTDist =
821  new TH1I("TofDigiBdf_ToTDist",
822  "ToT distribution, summed for all digis/channels; Tot [ns]",
823  500,
824  0.0,
825  50.0);
826 
827  fhElecChOccup = new TH1I("TofDigiBdf_ElecChOccup",
828  "Percentage of ToF electronic channels fired per "
829  "event; Elect. chan occupancy [%]",
830  1000,
831  0.0,
832  100.0);
833  fhMultiDigiEvtElCh =
834  new TH1D("TofDigiBdf_MultiDigiEvtElCh",
835  "Number of events with multiple digi (~MC track) per electronic "
836  "channel; Elect. chan index []",
837  fiNbElecChTot,
838  0.0,
839  fiNbElecChTot);
840  fhNbDigiEvtElCh =
841  new TH2D("TofDigiBdf_NbDigiEvtElCh",
842  "Number of digis per event before merging for each electronic "
843  "channel; Elect. chan index []; Nb Digis in Event []",
844  fiNbElecChTot,
845  0.0,
846  fiNbElecChTot,
847  10,
848  0.5,
849  10.5);
850  fhNbTracksEvtElCh =
851  new TH2D("TofDigiBdf_NbTracksEvtElCh",
852  "Number of tracks per event before merging for each electronic "
853  "channel; Elect. chan index []; Nb tracks in Event []",
854  fiNbElecChTot,
855  0.0,
856  fiNbElecChTot,
857  10,
858  0.5,
859  10.5);
860  fhFiredEvtElCh = new TH1D("TofDigiBdf_FiredEvtElCh",
861  "Number of events with at least 1 digi per "
862  "electronic channel; Elect. chan index []",
863  fiNbElecChTot,
864  0.0,
865  fiNbElecChTot);
866  fhMultiProbElCh = new TH1D(
867  "TofDigiBdf_MultiProbElCh",
868  "Probability of having multiple digi (~MC track) event per electronic "
869  "channel; Elect. chan index []; Multiple signal prob. [%]",
870  fiNbElecChTot,
871  0.0,
872  fiNbElecChTot);
873 
874  gDirectory->cd(
875  oldir
876  ->GetPath()); // <= To prevent histos from being sucked in by the param file of the TRootManager!
877 
878  return kTRUE;
879 }
880 Bool_t CbmTofDigitize::FillHistos() {
881 
882  // (VF) This does not run with digi vectors. Should be moved into a QA class.
883  /*
884  if(!fbMonitorHistos)
885  {
886  return kTRUE;
887  }
888 
889  Int_t nTofPoint = fTofPointsColl->GetEntries();
890  Int_t nTracks = fMcTracksColl->GetEntries();
891  // Int_t nTofDigi = fDigis->GetEntries();
892  Int_t iNbTofTracks = 0;
893  Double_t nTofFired = 0;
894  Double_t dProcessTime = fdDigitizeTime + fdMergeTime;
895 
896  // --- Update run Counters
897  fdNofPointsTot += nTofPoint;
898 
899  // Tracks Info
900  CbmMCTrack *pMcTrk;
901  if( fbMcTrkMonitor )
902  {
903  Int_t iNbTofTracksPrim = 0;
904  for(Int_t iTrkInd = 0; iTrkInd < nTracks; iTrkInd++)
905  {
906  pMcTrk = (CbmMCTrack*) fMcTracksColl->At( iTrkInd );
907  // Jump all tracks not making 8 points for test
908  // if( 8 != pMcTrk->GetNPoints(ECbmModuleId::kTof) )
909  // continue;
910  if( 0 < pMcTrk->GetNPoints(ECbmModuleId::kTof) )
911  {
912  iNbTofTracks++;
913  fhTofPointsPerTrack->Fill( pMcTrk->GetNPoints(ECbmModuleId::kTof) );
914  fhTofPointsPerTrackVsPdg->Fill( pMcTrk->GetPdgCode(), pMcTrk->GetNPoints(ECbmModuleId::kTof) );
915  }
916  if( 0 < pMcTrk->GetNPoints(ECbmModuleId::kTof) && -1 == pMcTrk->GetMotherId() )
917  iNbTofTracksPrim++;
918 
919  } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++)
920 
921  // --- Update run Counters
922  fdNofTofMcTrkTot += iNbTofTracks;
923  } // if( fbMcTrkMonitor )
924  LOG(info) << "Monitor tracks" ;
925 
926  // Tof Points info
927  CbmTofPoint *pTofPt;
928  for(Int_t iPtInd = 0; iPtInd < nTofPoint; iPtInd++)
929  {
930  pTofPt = (CbmTofPoint*) fTofPointsColl->At( iPtInd );
931  fhPtTime->Fill( pTofPt->GetTime() );
932 
933  if( fbMcTrkMonitor )
934  {
935  Int_t iDetId = pTofPt->GetDetectorID();
936  Int_t iGap = fGeoHandler->GetGap(iDetId);
937 
938  pMcTrk = (CbmMCTrack*) fMcTracksColl->At( pTofPt->GetTrackID() );
939 
940  fhTofPtsInTrkVsGapInd->Fill( pMcTrk->GetNPoints(ECbmModuleId::kTof), iGap );
941  if( -1 == pMcTrk->GetMotherId() )
942  fhTofPtsInTrkVsGapIndPrm->Fill( pMcTrk->GetNPoints(ECbmModuleId::kTof), iGap );
943  else if( 11 != pMcTrk->GetPdgCode() ) // Remove electrons
944  fhTofPtsInTrkVsGapIndSec->Fill( pMcTrk->GetNPoints(ECbmModuleId::kTof), iGap );
945 
946  // Get TofPoint Position
947  TVector3 vPntPos;
948  pTofPt->Position( vPntPos );
949  if( 8-pMcTrk->GetNPoints(ECbmModuleId::kTof) <= iGap &&
950  pMcTrk->GetNPoints(ECbmModuleId::kTof) < 8 &&
951  -1 != pMcTrk->GetMotherId() )
952  fhTofPtsPosVsGap[iGap]->Fill( vPntPos.X(), vPntPos.Y() );
953  } // if( fbMcTrkMonitor )
954  } // for(Int_t iPtInd = 0; iPtInd < nTofPoint; iPtInd++)
955  LOG(info) << "Monitor points" ;
956 
957 
958  fhDigiMergeTime->Fill( fdMergeTime );
959  fhEvtProcTime->Fill( dProcessTime );
960  fhProcTimeEvtSize->Fill( nTofPoint, dProcessTime );
961  if( fbMcTrkMonitor )
962  // fhMeanDigiPerTrack->Fill( (Double_t)nTofDigi/(Double_t)iNbTofTracks );
963 // fhMeanDigiPerPoint->Fill( (Double_t)nTofDigi/(Double_t)nTofPoint );
964 // fhMeanPointPerTrack->Fill( (Double_t)nTofPoint/(Double_t)iNbTofTracks);
965 // fhMeanPointPerTrack2d->Fill( iNbTofTracks, (Double_t)nTofPoint/(Double_t)iNbTofTracks);
966 
967  // Assume for the occupancy that we can only have one Digi per electronic channel
968  // No Multiple hits/digi in same event!
969  fhElecChOccup->Fill( 100.0*(Double_t)nTofDigi/(Double_t)fiNbElecChTot );
970  LOG(info) << "Monitor 1" ;
971 
972  if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
973  {
974  LOG(info) << "Monitor 2" ;
975  CbmTofDigiExp *pDigi;
976  for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
977  {
978  LOG(info) << "Digi " << iDigInd ;
979  pDigi = (CbmTofDigiExp*) fTofDigisColl->At( iDigInd );
980  assert(pDigi);
981  CbmMatch* digiMatch=(CbmMatch *)fTofDigiMatchPointsColl->At(iDigInd);
982  assert(digiMatch);
983  CbmLink L0 = digiMatch->GetMatchedLink();
984  CbmTofPoint* point = (CbmTofPoint*)fTofPointsColl->At( L0.GetIndex());
985  assert(point);
986 
987  if( pDigi->GetTot() < 0 )
988  cout<<iDigInd<<"/"<<nTofDigi<<" "<<pDigi->GetTot()<<endl;
989  fhDigiTime->Fill( pDigi->GetTime() );
990  fhDigiTimeRes->Fill( pDigi->GetTime() - point->GetTime() );
991  fhDigiTimeResB->Fill( pDigi->GetTime() - point->GetTime(),
992  pDigi->GetType() );
993  fhToTDist->Fill( pDigi->GetTot() );
994 
995  nTofFired += 1.0/( 2.0 - fDigiBdfPar->GetChanType( pDigi->GetType(), pDigi->GetRpc() ) );
996  } // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
997  } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
998  else
999  {
1000  LOG(info) << "Monitor 3" ;
1001 
1002  CbmTofDigi *pDigi;
1003  for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
1004  {
1005  pDigi = (CbmTofDigi*) fTofDigisColl->At( iDigInd );
1006  CbmMatch* digiMatch=(CbmMatch *)fTofDigiMatchPointsColl->At(iDigInd);
1007  CbmLink L0 = digiMatch->GetMatchedLink();
1008  CbmTofPoint* point = (CbmTofPoint*)fTofPointsColl->At( L0.GetIndex());
1009  // CbmTofPoint* point = (CbmTofPoint*)fTofPointsColl->At(pDigi->GetMatch()->GetMatchedLink().GetIndex());
1010  if( pDigi->GetTot() < 0 )
1011  cout<<iDigInd<<"/"<<nTofDigi<<" "<<pDigi->GetTot()<<endl;
1012  fhDigiTime->Fill( pDigi->GetTime() );
1013  fhDigiTimeRes->Fill( pDigi->GetTime() - point->GetTime() );
1014  fhToTDist->Fill( pDigi->GetTot() );
1015 
1016  nTofFired += 1.0/( 2.0 - fDigiBdfPar->GetChanType( pDigi->GetType(), pDigi->GetRpc() ) );
1017  } // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
1018  } // else of if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
1019 
1020 // fhMeanFiredPerPoint->Fill( nTofFired / (Double_t)nTofPoint);
1021  fhMeanFiredPerTrack->Fill( nTofFired/(Double_t)iNbTofTracks );
1022  LOG(info) << "Monitor time" ;
1023 */
1024 
1025  return kTRUE;
1026 }
1027 Bool_t CbmTofDigitize::SetHistoFileName(TString sFilenameIn) {
1028  fsHistoOutFilename = sFilenameIn;
1029  return kTRUE;
1030 }
1031 
1032 
1033 Bool_t CbmTofDigitize::WriteHistos() {
1034  if ("" == fsHistoOutFilename || !fbMonitorHistos) {
1035  LOG(info) << "CbmTofDigitize::WriteHistos => Control histograms will not "
1036  "be written to disk!";
1037  LOG(info) << "CbmTofDigitize::WriteHistos => To change this behavior "
1038  "please provide a full path "
1039  << "with the SetHistoFileName method";
1040  return kTRUE;
1041  } // if( "" == fsHistoOutFilename )
1042 
1043  TDirectory* oldir = gDirectory;
1044  TFile* fHist = new TFile(fsHistoOutFilename, "RECREATE");
1045 
1046  fHist->cd();
1047 
1048  fhTofPointsPerTrack->Write();
1049  fhTofPtsInTrkVsGapInd->Write();
1050  fhTofPtsInTrkVsGapIndPrm->Write();
1051  fhTofPtsInTrkVsGapIndSec->Write();
1052  for (Int_t iGap = 0; iGap < 10; iGap++)
1053  fhTofPtsPosVsGap[iGap]->Write();
1054  /*
1055  fhTofPointsPerTrackVsPdg->Write();
1056  fhMeanDigiPerPoint->Write();
1057  fhMeanFiredPerPoint->Write();
1058  fhMeanPointPerTrack->Write();
1059  fhMeanPointPerTrack2d->Write();
1060 */
1061  fhEvtProcTime->Write();
1062  fhDigiMergeTime->Write();
1063  fhDigiNbElecCh->Write();
1064  fhProcTimeEvtSize->Write();
1065  fhMeanDigiPerTrack->Write();
1066  fhMeanFiredPerTrack->Write();
1067  fhPtTime->Write();
1068  fhDigiTime->Write();
1069  fhDigiTimeRes->Write();
1070  fhDigiTimeResB->Write();
1071  fhToTDist->Write();
1072 
1073  fhElecChOccup->Write();
1074  fhMultiDigiEvtElCh->Write();
1075  fhNbDigiEvtElCh->Write();
1076  fhNbTracksEvtElCh->Write();
1077  fhFiredEvtElCh->Write();
1078  fhMultiProbElCh->Divide(fhMultiDigiEvtElCh, fhFiredEvtElCh);
1079  fhMultiProbElCh->Scale(100.0);
1080  fhMultiProbElCh->Write();
1081 
1082  // Uncomment if need to control the Cluster Size and ToT proba
1083  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
1084  for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
1085  if (0 == fDigiBdfPar->GetClusterModel())
1086  fh1ClusterSizeProb[iSmType]->Write();
1087  fh1ClusterTotProb[iSmType]->Write();
1088  }
1089 
1090  gDirectory->cd(oldir->GetPath());
1091 
1092  fHist->Close();
1093 
1094  return kTRUE;
1095 }
1096 Bool_t CbmTofDigitize::DeleteHistos() {
1097  if (!fbMonitorHistos) { return kTRUE; }
1098 
1099  delete fhTofPointsPerTrack;
1100  delete fhTofPtsInTrkVsGapInd;
1101  delete fhTofPtsInTrkVsGapIndPrm;
1102  delete fhTofPtsInTrkVsGapIndSec;
1103  for (Int_t iGap = 0; iGap < 10; iGap++)
1104  delete fhTofPtsPosVsGap[iGap];
1105  /*
1106  delete fhTofPointsPerTrackVsPdg;
1107  delete fhMeanDigiPerPoint;
1108  delete fhMeanFiredPerPoint;
1109  delete fhMeanPointPerTrack;
1110  delete fhMeanPointPerTrack2d;
1111 */
1112  delete fhEvtProcTime;
1113  delete fhDigiMergeTime;
1114  delete fhDigiNbElecCh;
1115  delete fhProcTimeEvtSize;
1116  delete fhMeanDigiPerTrack;
1117  delete fhMeanFiredPerTrack;
1118  delete fhPtTime;
1119  delete fhDigiTime;
1120  delete fhDigiTimeRes;
1121  delete fhDigiTimeResB;
1122  delete fhToTDist;
1123 
1124  delete fhElecChOccup;
1125  delete fhMultiDigiEvtElCh;
1126  delete fhNbDigiEvtElCh;
1127  delete fhNbTracksEvtElCh;
1128  delete fhFiredEvtElCh;
1129  delete fhMultiProbElCh;
1130 
1131  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
1132  for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
1133  if (0 == fDigiBdfPar->GetClusterModel()) delete fh1ClusterSizeProb[iSmType];
1134  delete fh1ClusterTotProb[iSmType];
1135  }
1136 
1137  return kTRUE;
1138 }
1139 /************************************************************************************/
1140 // Functions for the merging of "gap digis" and "multiple hits digis" into "channel digis"
1141 // TODO: FEE double hit discrimination capability (using Time distance between Digis)
1142 // TODO: Charge summing up
1143 Bool_t CbmTofDigitize::MergeSameChanDigis() {
1144 
1145  if (bFakeBeamCounter) {
1146  UInt_t uChanUId = 0x00005006;
1147  Double_t dHitTime = fCurrentEventTime + gRandom->Gaus(0., 0.04);
1148  const Double_t dHitTot = 2.;
1149  CbmTofDigi* tDigi = new CbmTofDigi(uChanUId, dHitTime, dHitTot);
1150  CbmMatch* tMatch = new CbmMatch();
1151  SendData(tDigi, tMatch); // Send digi to DAQ
1152  fiNbDigis++;
1153  LOG(debug) << Form(
1154  "Add fake diamond digis 0x%08x mode with t = %7.3f", uChanUId, dHitTime);
1155  //delete tDigi;
1156  //delete tMatch;
1157  }
1158 
1159  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
1160 
1161  // loop over each (Smtype, Sm, Rpc, Channel, Side)
1162  for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
1163  Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
1164  Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
1165 
1166  for (Int_t iSm = 0; iSm < iNbSm; iSm++)
1167  for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
1168  Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
1169  Int_t iNbSides = 2 - fDigiBdfPar->GetChanType(iSmType, iRpc);
1170 
1171  for (Int_t iCh = 0; iCh < iNbCh; iCh++)
1172  for (Int_t iSide = 0; iSide < iNbSides; iSide++) {
1173  Int_t iNbDigis =
1174  fStorDigi[iSmType][iSm * iNbRpc + iRpc][iNbSides * iCh + iSide]
1175  .size();
1176 
1177  // TOF QA: monitor number of tracks leading to digi before merging them
1178  Int_t iNbTracks = 0;
1179  std::vector<Int_t> vPrevTrackIdList;
1180  if (0 < iNbDigis) {
1181  CbmMatch* digiMatch;
1182  for (Int_t iDigi = 0; iDigi < iNbDigis;
1183  iDigi++) { //store matches with digis
1184  LOG(DEBUG) << Form(
1185  " Create match for Digi %d(%d), addr 0x%08x; %d",
1186  iDigi,
1187  iNbDigis,
1188  fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1189  [iNbSides * iCh + iSide][iDigi]
1190  .first->GetAddress(),
1191  fStorDigiMatch[iSmType][iSm * iNbRpc + iRpc]
1192  [iNbSides * iCh + iSide][iDigi]);
1193 
1194  digiMatch = new CbmMatch();
1195  digiMatch->AddLink(
1196  CbmLink(1.,
1197  fStorDigiMatch[iSmType][iSm * iNbRpc + iRpc]
1198  [iNbSides * iCh + iSide][iDigi],
1199  fCurrentMCEntry,
1200  fCurrentInput));
1201 
1202  fStorDigi[iSmType][iSm * iNbRpc + iRpc][iNbSides * iCh + iSide]
1203  [iDigi]
1204  .second = digiMatch;
1205  }
1206  if (1 < iNbDigis) { // time sort digis
1207  std::sort(fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1208  [iNbSides * iCh + iSide]
1209  .begin(),
1210  fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1211  [iNbSides * iCh + iSide]
1212  .end(),
1213  CompTExp);
1214  }
1215 
1216  if (fbMonitorHistos) {
1217  fhNbDigiEvtElCh->Fill(fvRpcChOffs[iSmType][iSm][iRpc]
1218  + iNbSides * iCh + iSide,
1219  iNbDigis);
1220 
1221  fhDigiNbElecCh->Fill(iNbDigis);
1222 
1223  fhFiredEvtElCh->Fill(fvRpcChOffs[iSmType][iSm][iRpc]
1224  + iNbSides * iCh + iSide);
1225  if (1 < iNbDigis)
1226  fhMultiDigiEvtElCh->Fill(fvRpcChOffs[iSmType][iSm][iRpc]
1227  + iNbSides * iCh + iSide);
1228  }
1229 
1230  if (0
1231  == fStorDigiMatch[iSmType][iSm * iNbRpc + iRpc]
1232  [iNbSides * iCh + iSide]
1233  .size()) {
1234  LOG(error) << Form(
1235  " cannot add digiMatch for (%d,%d,%d,%d,%d) at pos %d",
1236  iSmType,
1237  iSm,
1238  iRpc,
1239  iCh,
1240  iSide,
1241  fiNbDigis);
1242  break;
1243  }
1244 
1245  Int_t iDigi0 = 0;
1246  digiMatch = fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1247  [iNbSides * iCh + iSide][iDigi0]
1248  .second;
1249  if (iNbDigis > 1) {
1250  LOG(debug) << Form("Now apply deadtime %f for %d digis, "
1251  "digi0=%d for adress 0x%08x ",
1252  fDigiBdfPar->GetDeadtime(),
1253  iNbDigis,
1254  iDigi0,
1255  fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1256  [iNbSides * iCh + iSide][iDigi0]
1257  .first->GetAddress());
1258  Double_t fdDeadtime = fDigiBdfPar->GetDeadtime();
1259  for (Int_t iDigi = 1; iDigi < iNbDigis; iDigi++) {
1260  // find Digis in deadtime of first Digi
1261  if (
1262  (fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1263  [iNbSides * iCh + iSide][iDigi]
1264  .first->GetTime()
1265  - fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1266  [iNbSides * iCh + iSide][iDigi0]
1267  .first->GetTime())
1268  < fdDeadtime) { // Store Link with weight 0 to have a trace of MC tracks firing the channel
1269  // but whose Digi is not propagated to next stage
1270  digiMatch = fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1271  [iNbSides * iCh + iSide][iDigi0]
1272  .second;
1273  Int_t nL = fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1274  [iNbSides * iCh + iSide][iDigi]
1275  .second->GetNofLinks();
1276  if (nL != 1) LOG(fatal) << "Invalid number of Links " << nL;
1277 
1278  CbmLink L0 = fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1279  [iNbSides * iCh + iSide][iDigi]
1280  .second->GetLink(0);
1281  digiMatch->AddLink(CbmLink(
1282  0., L0.GetIndex(), fCurrentMCEntry, fCurrentInput));
1283 
1284  // cleanup of obsolete digis, matches and links done below
1285  } else { // new digi found
1286  // new digi will be created and later deleted by CbmDaq.
1287  // The original digi will be deleted below, together with the unused digis from the buffer.
1288  CbmTofDigi* digi =
1289  new CbmTofDigi(*(fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1290  [iNbSides * iCh + iSide][iDigi0]
1291  .first));
1292  CbmMatch* match =
1293  new CbmMatch(*(fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1294  [iNbSides * iCh + iSide][iDigi0]
1295  .second));
1296 
1297  digi->SetTime(digi->GetTime() * fdDigiTimeConvFactor
1298  + fCurrentEventTime); // ns->ps
1299  SendData(digi, match); // Send digi to DAQ
1300  fiNbDigis++;
1301 
1302  // TOF QA
1303  if (fbMonitorHistos && NULL != digiMatch) {
1304  Int_t nL = digiMatch->GetNofLinks();
1305  for (Int_t iL = 0; iL < nL; iL++) {
1306  // TOF QA: monitor number of tracks leading to digi before merging them
1307  CbmLink L = digiMatch->GetLink(iL);
1308  Bool_t bTrackFound = kFALSE;
1309  for (UInt_t uTrkId = 0;
1310  uTrkId < vPrevTrackIdList.size();
1311  uTrkId++)
1312  if (vPrevTrackIdList[uTrkId]
1313  == ((CbmTofPoint*) fTofPointsColl->At(
1314  L.GetIndex()))
1315  ->GetTrackID()) {
1316  bTrackFound = kTRUE;
1317  break;
1318  } // If TrackId already found for this Elec. Chan in this event
1319  if (kFALSE == bTrackFound) {
1320  // New track or first track
1321  iNbTracks++;
1322  vPrevTrackIdList.push_back(
1323  ((CbmTofPoint*) fTofPointsColl->At(L.GetIndex()))
1324  ->GetTrackID());
1325  } // if( kFALSE == bTrackFound )
1326  } // for( Int_t iL = 0; iL < nL; iL++
1327  vPrevTrackIdList.clear();
1328  // TOF QA: monitor number of tracks leading to digi before merging them
1329  fhNbTracksEvtElCh->Fill(fvRpcChOffs[iSmType][iSm][iRpc]
1330  + iNbSides * iCh + iSide,
1331  iNbTracks);
1332  } // if(fbMonitorHistos)
1333 
1334  iDigi0 = iDigi;
1335  }
1336  }
1337  } // if( iNbDigis > 1) end
1338  // assemble and send last digi
1339  // new digi will be created and later deleted by CbmDaq.
1340  // The original digi will be deleted below, together with the unused digis from the buffer.
1341  CbmTofDigi* digi =
1342  new CbmTofDigi(*(fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1343  [iNbSides * iCh + iSide][iDigi0]
1344  .first));
1345  CbmMatch* match =
1346  new CbmMatch(*(fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1347  [iNbSides * iCh + iSide][iDigi0]
1348  .second));
1349  digi->SetTime(digi->GetTime() * fdDigiTimeConvFactor
1350  + fCurrentEventTime); // ns->ps
1351  SendData(digi, match); // Send digi to DAQ
1352  fiNbDigis++;
1353 
1354  if (fbMonitorHistos) {
1355  Int_t nL = digiMatch->GetNofLinks();
1356  for (Int_t iL = 0; iL < nL; iL++) {
1357  // TOF QA: monitor number of tracks leading to digi before merging them
1358  CbmLink L = digiMatch->GetLink(iL);
1359  Bool_t bTrackFound = kFALSE;
1360  for (UInt_t uTrkId = 0; uTrkId < vPrevTrackIdList.size();
1361  uTrkId++)
1362  if (vPrevTrackIdList[uTrkId]
1363  == ((CbmTofPoint*) fTofPointsColl->At(L.GetIndex()))
1364  ->GetTrackID()) {
1365  bTrackFound = kTRUE;
1366  break;
1367  } // If TrackId already found for this Elec. Chan in this event
1368  if (kFALSE == bTrackFound) {
1369  // New track or first track
1370  iNbTracks++;
1371  vPrevTrackIdList.push_back(
1372  ((CbmTofPoint*) fTofPointsColl->At(L.GetIndex()))
1373  ->GetTrackID());
1374  } // if( kFALSE == bTrackFound )
1375  } // for( Int_t iL = 0; iL < nL; iL++
1376  vPrevTrackIdList.clear();
1377  // TOF QA: monitor number of tracks leading to digi before merging them
1378  fhNbTracksEvtElCh->Fill(fvRpcChOffs[iSmType][iSm][iRpc]
1379  + iNbSides * iCh + iSide,
1380  iNbTracks);
1381  } // if(fbMonitorHistos)
1382 
1383  for (Int_t iDigi = 0; iDigi < iNbDigis; iDigi++) {
1384  CbmTofDigi* tmpDigi = fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1385  [iNbSides * iCh + iSide][iDigi]
1386  .first;
1387  CbmMatch* tmpMatch = fStorDigi[iSmType][iSm * iNbRpc + iRpc]
1388  [iNbSides * iCh + iSide][iDigi]
1389  .second;
1390  if (tmpDigi) delete tmpDigi;
1391  if (tmpMatch) delete tmpMatch;
1392  // delete fStorDigiMatch[iSmType][iSm*iNbRpc + iRpc][iNbSides*iCh+iSide][iDigi];
1393  } //
1394 
1395  fStorDigi[iSmType][iSm * iNbRpc + iRpc][iNbSides * iCh + iSide]
1396  .clear();
1397  fStorDigiMatch[iSmType][iSm * iNbRpc + iRpc]
1398  [iNbSides * iCh + iSide]
1399  .clear();
1400  } // if( 0 < iNbDigis )
1401  } // for each (Ch, Side) pair
1402  } // for each (Sm, rpc) pair
1403  } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
1404 
1405  return kTRUE;
1406 }
1407 /************************************************************************************/
1408 // Functions for the Cluster Radius generation
1409 Double_t CbmTofDigitize::GenerateClusterRadius(Int_t iSmType, Int_t iRpc) {
1410  Double_t dClusterRadius = 0;
1411  Int_t iChType = fDigiBdfPar->GetChanType(iSmType, iRpc);
1412 
1413  switch (fDigiBdfPar->GetClusterRadiusModel()) {
1414  case -1: {
1415  // Fixed value at 0.0002 to get a cluster size as close to 1 as possible
1416  dClusterRadius = 0.0002;
1417  break;
1418  }
1419  case 0: {
1420  // Single value using the beamtime cluster size distribution mean
1421  if (0 == iChType) {
1422  // Simple linear relation mean cluster size = 2* radius +1 in [strips]
1423  // Cf "RadiusToMeanClusterAll" histogram
1424  // Come from simple geometry of cluster center position if one
1425  // neglect RPC and extremities edge effects
1426  dClusterRadius = (fh1ClusterSizeProb[iSmType]->GetMean() - 1.0) / 2.0;
1427  // Convert to [cm]
1428  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
1429  // Horizontal channel => strip width along Y
1430  dClusterRadius *= fChannelInfo->GetSizey();
1431  // Vertical channel => strip width along X
1432  else
1433  dClusterRadius *= fChannelInfo->GetSizex();
1434  } // if( 0 == iChType)
1435  else {
1436  LOG(error)
1437  << "CbmTofDigitize::GenerateClusterRadius => Cluster Radius "
1438  << "obtention from cluster size not implemented for pads, Sm type "
1439  << iSmType << " Rpc " << iRpc;
1440  LOG(error) << "CbmTofDigitize::GenerateClusterRadius => Test "
1441  << " Sm type " << iSmType << " Rpc " << iRpc << " Type "
1442  << iChType << " Type " << (Int_t) iChType;
1443  return kFALSE;
1444  } // else of if( 0 == iChType)
1445  break;
1446  } // case 0:
1447  case 1:
1448  case 2: {
1449  if (0 == iChType) {
1450  // from beamtime cluster size distribution
1451  Double_t dStripSize = 0;
1452  // Convert to [cm]
1453  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
1454  // Horizontal channel => strip width along Y
1455  dStripSize = fChannelInfo->GetSizey();
1456  // Vertical channel => strip width along X
1457  else
1458  dStripSize = fChannelInfo->GetSizex();
1459 
1460  // Obtain the Landau peak position and sigma scale factor
1461  // from the parameter object. It should take care of getting the values
1462  // either from the ASCII parameter file or from a fit on the beam data
1463  // TODO: cross-check the exact value of the second parameter, here it is set so that probability of
1464  // negative radius is really low
1465  Double_t dPeakPos = fDigiBdfPar->GetLandauMpv(iSmType);
1466  Double_t dSigmScal =
1467  fDigiBdfPar->GetLandauSigma(iSmType); // empirical best
1468 
1469  // Get the actual Cluster radius
1470  dClusterRadius = dStripSize * gRandom->Landau(dPeakPos, dSigmScal);
1471  if (dClusterRadius < 0) dClusterRadius = 0;
1472  } // if( 0 == iChType)
1473  else {
1474  LOG(error)
1475  << "CbmTofDigitize::GenerateClusterRadius => Cluster Radius "
1476  << "obtention from cluster size not implemented for pads, Sm type "
1477  << iSmType << " Rpc " << iRpc;
1478  LOG(error) << "CbmTofDigitize::DigitizeFlatDisc => Test 2"
1479  << " Sm type " << iSmType << " Rpc " << iRpc << " Type "
1480  << iChType;
1481  return kFALSE;
1482  } // else of if( 0 == iChType)
1483  break;
1484  } // case 1:
1485  default: {
1486  LOG(error) << "CbmTofDigitize::GenerateClusterRadius => Wrong Cluster "
1487  "Radius method , Sm type "
1488  << iSmType << " Rpc " << iRpc;
1489  dClusterRadius = 0;
1490  break;
1491  } // default:
1492  } // switch( fDigiBdfPar->GetClusterRadiusModel() )
1493  return dClusterRadius;
1494 }
1495 /************************************************************************************/
1496 // Functions for a direct use of the cluster size
1497 Bool_t CbmTofDigitize::DigitizeDirectClusterSize() {
1498  // Uniform distribution in ]0;x]
1499  // gRandom->Uniform(x);
1500  // Gauss distribution in of mean m, sigma s
1501  // gRandom->Gauss(m, s);
1502 
1503  CbmTofPoint* pPoint;
1504  CbmMCTrack* pMcTrk;
1505 
1506  Int_t nTofPoint = fTofPointsColl->GetEntries();
1507  Int_t nMcTracks = fMcTracksColl->GetEntries();
1508 
1509  // Prepare the temporary storing of the Track/Point/Digi info
1510  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
1511  fvlTrckChAddr.resize(nMcTracks);
1512  fvlTrckRpcAddr.resize(nMcTracks);
1513  fvlTrckRpcTime.resize(nMcTracks);
1514  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
1515  for (Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++) {
1516  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
1517  fvlTrckChAddr[iTrkInd].clear();
1518  fvlTrckRpcAddr[iTrkInd].clear();
1519  fvlTrckRpcTime[iTrkInd].clear();
1520  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
1521  } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++)
1522 
1523  if (fbMcTrkMonitor) {
1524  // Debug Info printout
1525  Int_t iNbTofTracks = 0;
1526  Int_t iNbTofTracksPrim = 0;
1527 
1528  LOG(debug1) << "CbmTofDigitize::DigitizeDirectClusterSize: " << nTofPoint
1529  << " points in Tof for this event with " << nMcTracks
1530  << " MC tracks ";
1531  for (Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++) {
1532  pMcTrk = (CbmMCTrack*) fMcTracksColl->At(iTrkInd);
1533  if (0 < pMcTrk->GetNPoints(ECbmModuleId::kTof)) iNbTofTracks++;
1534  if (0 < pMcTrk->GetNPoints(ECbmModuleId::kTof)
1535  && -1 == pMcTrk->GetMotherId())
1536  iNbTofTracksPrim++;
1537  } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++)
1538 
1539  //Some numbers on TOF distributions
1540  LOG(debug1) << "CbmTofDigitize::DigitizeDirectClusterSize: " << iNbTofTracks
1541  << " tracks in Tof ";
1542  LOG(debug1) << "CbmTofDigitize::DigitizeDirectClusterSize: "
1543  << iNbTofTracksPrim << " tracks in Tof from vertex";
1544  } // if( fbMcTrkMonitor )
1545 
1546  // Tof Point Info
1547  Int_t iDetId, iSmType, iSM, iRpc, iChannel, iGap;
1548  Int_t iTrackID, iChanId;
1549  TVector3 vPntPos;
1550  // Cluster Info
1551  Int_t iClusterSize;
1552  Double_t dClustCharge;
1553  // Digi
1554  // CbmTofDigi * pTofDigi; // <- Compressed digi (64 bits)
1555  // CbmTofDigi * pTofDigiExp; // <- Expanded digi
1556  // General geometry info
1557  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
1558  Int_t iNbSm, iNbRpc, iNbCh;
1559  Int_t iChType;
1560 
1561  // Start jitter: Same contrib. for all points in same event
1562  // FIXME
1563  // Actually, event start time reconstruction should not be done in the ToF
1564  // digitizer, neither in event-based nor in time-based mode. The timestamp
1565  // of a CbmTofDigi object is not the time difference between signals in a
1566  // (hypothetical) start detector system and the MRPC wall but rather
1567  // represents the detection point in time in an MRPC measured with respect to
1568  // the start of the run.
1569  // For compatibility reasons, we continue adding a start detector jitter to
1570  // digi timestamps in event-based mode. In time-based mode, event start time
1571  // reconstruction is not considered here.
1572  Double_t dStartJitter = 0.;
1573  if (fEventMode) {
1574  dStartJitter = gRandom->Gaus(0.0, fDigiBdfPar->GetStartTimeRes());
1575  }
1576 
1577  for (Int_t iPntInd = 0; iPntInd < nTofPoint; iPntInd++) {
1578  // Get a pointer to the TOF point
1579  pPoint = (CbmTofPoint*) fTofPointsColl->At(iPntInd);
1580  if (NULL == pPoint) {
1581  LOG(WARNING) << "CbmTofDigitize::DigitizeDirectClusterSize => Be "
1582  "careful: hole in the CbmTofPoint TClonesArray!";
1583  continue;
1584  } // if( pPoint )
1585 
1586  // Get all channel info
1587  iDetId = pPoint->GetDetectorID();
1588  iSmType = fGeoHandler->GetSMType(iDetId);
1589  iRpc = fGeoHandler->GetCounter(iDetId);
1590 
1591  iChannel = fGeoHandler->GetCell(iDetId);
1592  if (fGeoHandler->GetGeoVersion() < k14a)
1593  iChannel--; // Again, channel going from 1 to nbCh instead of 0 to nbCh - 1 ?!?!?
1594  iGap = fGeoHandler->GetGap(iDetId);
1595  iSM = fGeoHandler->GetSModule(iDetId);
1596  iChanId = fGeoHandler->GetCellId(iDetId);
1597  iTrackID = pPoint->GetTrackID();
1598 
1599  iNbSm = fDigiBdfPar->GetNbSm(iSmType);
1600  iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
1601  iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
1602  iChType = fDigiBdfPar->GetChanType(iSmType, iRpc);
1603  Int_t iTrkId = pPoint->GetTrackID();
1604 
1605  // Catch case where the MC track was removed from the transport stack
1606  if (iTrkId < 0) {
1607  if (fbAllowPointsWithoutTrack)
1608  continue;
1609  else
1610  LOG(fatal)
1611  << "CbmTofDigitize::DigitizeDirectClusterSize => TofPoint without "
1612  "valid MC track Index, "
1613  << " track was probably cut at the transport level! "
1614  << " Pnt Idx: " << iPntInd << " Trk Idx: " << iTrkId << "\n"
1615  << "=============> To allow this kind of points and simply jump "
1616  "them, "
1617  << "call CbmTofDigitize::AllowPointsWithoutTrack() in your macro!!";
1618  } // if( iTrkId < 0 )
1619 
1620  if (fbMcTrkMonitor) {
1621  // Get pointer to the MC-Track info
1622  pMcTrk = (CbmMCTrack*) fMcTracksColl->At(iTrkId);
1623 
1624  // Keep only primaries
1625  if (kTRUE == fDigiBdfPar->UseOnlyPrimaries())
1626  if (-1 != pMcTrk->GetMotherId()) continue;
1627  } // if( fbMcTrkMonitor )
1628 
1629  if (iSmType < 0. || iNbSmTypes <= iSmType || iSM < 0. || iNbSm <= iSM
1630  || iRpc < 0. || iNbRpc <= iRpc || iChannel < 0. || iNbCh <= iChannel) {
1631  LOG(error) << Form("CbmTofDigitize => det ID 0x%08x", iDetId)
1632  << " SMType: " << iSmType;
1633  LOG(error) << " SModule: " << iSM << " of " << iNbSm + 1;
1634  LOG(error) << " Module: " << iRpc << " of " << iNbRpc + 1;
1635  LOG(error) << " Gap: " << iGap;
1636  LOG(fatal) << " Cell: " << iChannel << " of " << iNbCh + 1;
1637  continue;
1638  } // if error on channel data
1639 
1640  // Check if there was already a Digi from the same track created
1641  // with this channel as main channel
1642  ULong64_t uAddr =
1643  CbmTofAddress::GetUniqueAddress(iSM, iRpc, iChannel, 0, iSmType);
1644  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
1645  Bool_t bFoundIt = kFALSE;
1646  // Check if this track ID is bigger than size of McTracks TClonesArray (maybe happen with PSD cut!!)
1647  // If yes, resize the array to reach the appropriate number of fields
1648  // Cast of TrackId is safe as we already check for "< 0" case
1649  // TODO: Is this check still required when Id < 0 are jumped?
1650  if (fvlTrckChAddr.size() <= static_cast<UInt_t>(iTrkId))
1651  fvlTrckChAddr.resize(iTrkId + 1);
1652 
1653  for (UInt_t uTrkMainCh = 0; uTrkMainCh < fvlTrckChAddr[iTrkId].size();
1654  uTrkMainCh++)
1655  if (uAddr == fvlTrckChAddr[iTrkId][uTrkMainCh]) {
1656  bFoundIt = kTRUE;
1657  break;
1658  }
1659  // If it is the case, no need to redo the full stuff for this gap
1660  if (kTRUE == bFoundIt) continue;
1661  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
1662 
1663  // Probability that the RPC detect the hit
1664  // if( fDigiBdfPar->GetEfficiency(iSmType) < gRandom->Uniform(1) )
1665  // continue;
1666 
1667  // Probability that the gap detect the hit
1668  // For now use a simple gap treatment (cf CbmTofDigiBdfPar):
1669  // - a single fired gap fires the channel
1670  // - all gaps have exactly equivalent efficiency/firing probability
1671  // - independent gap firing (no charge accumulation or such)
1672  // The probability for a hit not to fire at all is then
1673  // (1-p)^NGaps with p = gap efficiency and Ngaps the number of gaps in this RPC
1674  // This results in the relation: p = 1 - (1 - P)^(1/Ngaps)
1675  // with P = RPC efficiency from beamtime
1676  if (fDigiBdfPar->GetGapEfficiency(iSmType, iRpc) < gRandom->Uniform(1))
1677  continue;
1678 
1679  // Save the address in vector so that this channel is not redone later for the
1680  // eventual other gaps TofPoint
1681  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk())
1682  fvlTrckChAddr[iTrkId].push_back(uAddr);
1683 
1684  // Get TofPoint Position
1685  pPoint->Position(vPntPos);
1686  fChannelInfo = fDigiPar->GetCell(iChanId);
1687  TGeoNode* fNode = // prepare global->local trafo
1688  gGeoManager->FindNode(
1689  fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
1690  LOG(debug1) << Form(" TofDigitizerBDF:: (%3d,%3d,%3d,%3d) - node at "
1691  "(%6.1f,%6.1f,%6.1f) : 0x%p Pos(%6.1f,%6.1f,%6.1f)",
1692  iSmType,
1693  iSM,
1694  iRpc,
1695  iChannel,
1696  fChannelInfo->GetX(),
1697  fChannelInfo->GetY(),
1698  fChannelInfo->GetZ(),
1699  fNode,
1700  vPntPos.X(),
1701  vPntPos.Y(),
1702  vPntPos.Z());
1703  gGeoManager->GetCurrentNode();
1704  gGeoManager->GetCurrentMatrix();
1705 
1706  Double_t poipos[3] = {vPntPos.X(), vPntPos.Y(), vPntPos.Z()};
1707  Double_t poipos_local[3];
1708  Double_t poipos_local_man[3];
1709  gGeoManager->MasterToLocal(poipos, poipos_local_man);
1710  fDigiPar->GetNode(iChanId)->MasterToLocal(poipos, poipos_local);
1711  for (Int_t i = 0; i < 3; i++)
1712  if (poipos_local[i] != poipos_local_man[i])
1713  LOG(fatal) << "Inconsistent M2L result " << i << ": " << poipos_local[i]
1714  << " != " << poipos_local_man[i];
1715 
1716  if (1 == iChType) {
1717  LOG(error) << "CbmTofDigitize::DigitizeDirectClusterSize => This method "
1718  << "is not available for pads!!";
1719  return kFALSE;
1720  } // if( 1 == iChType)
1721 
1722  // Cluster size from beamtime distribution in [strips]
1723  iClusterSize = fh1ClusterSizeProb[iSmType]->GetBinContent(
1724  fh1ClusterSizeProb[iSmType]->FindBin(gRandom->Uniform(1)));
1725 
1726  Int_t iNbStripsSideA;
1727  Int_t iNbStripsSideB;
1728  if (1 == iClusterSize % 2) {
1729  // odd => a center strip and symetric side strips
1730  iNbStripsSideA = (iClusterSize - 1) / 2;
1731  iNbStripsSideB = (iClusterSize - 1) / 2;
1732  } // if odd cluster size
1733  else {
1734  // even => asymetric, the side getting more strips depends
1735  // on the cluster center position
1736  iNbStripsSideA = iClusterSize / 2; // left/bottom
1737  iNbStripsSideB = iClusterSize / 2; // right/top
1738 
1739  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
1740  // Horizontal strips
1741  if (poipos_local[1]
1742  < fChannelInfo->GetSizey()
1743  / 2.0) //FIXME, needs posiion in rotated counter frame
1744  iNbStripsSideB--; // less strips on top
1745  else
1746  iNbStripsSideA--; // less strips on bottom
1747  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
1748  else if (poipos_local[0] < fChannelInfo->GetSizex() / 2.0) //FIXME
1749  iNbStripsSideB--; // less strips on right
1750  else
1751  iNbStripsSideA--; // less strips on left
1752  } // if even cluster size
1753 
1754  // Generate Cluster charge as ToT[ps]
1755  dClustCharge = fh1ClusterTotProb[iSmType]->GetBinContent(
1756  fh1ClusterTotProb[iSmType]->FindBin(gRandom->Uniform(1)));
1757 
1758  // Calculate the time for the central channel
1759  Double_t dCentralTime;
1760  // Check if there was already a Digi from the same track created in this RPC
1761  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
1762  ULong64_t uRpcAddr =
1763  CbmTofAddress::GetUniqueAddress(iSM, iRpc, 0, 0, iSmType);
1764  Bool_t bFoundIt = kFALSE;
1765  UInt_t uTrkRpcPair = 0;
1766  for (uTrkRpcPair = 0; uTrkRpcPair < fvlTrckRpcAddr[iTrkId].size();
1767  uTrkRpcPair++)
1768  if (uRpcAddr == fvlTrckRpcAddr[iTrkId][uTrkRpcPair]) {
1769  bFoundIt = kTRUE;
1770  break;
1771  }
1772  // If it is the case, we should reuse the timing already assigned to this track
1773  if (kTRUE == bFoundIt)
1774  dCentralTime = fvlTrckRpcTime[iTrkId][uTrkRpcPair];
1775  else {
1776  dCentralTime =
1777  pPoint->GetTime()
1778  + gRandom->Gaus(0.0, fDigiBdfPar->GetResolution(iSmType))
1779  + dStartJitter; // Same contrib. for all points in same event
1780  fvlTrckRpcAddr[iTrkId].push_back(uRpcAddr);
1781  fvlTrckRpcTime[iTrkId].push_back(dCentralTime);
1782  } // No Digi yet in this RPC for this Track
1783  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
1784  else
1785  dCentralTime =
1786  pPoint->GetTime()
1787  + gRandom->Gaus(0.0, fDigiBdfPar->GetResolution(iSmType))
1788  + dStartJitter; // Same contrib. for all points in same event
1789 
1790  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
1791  // Horizontal channel: A = Right and B = Left of the channel
1792 
1793  // Assume the width (sigma) of the gaussian is 0.5*dClusterSize/3
1794  // => 99% of the charge is in "dClusterSize"
1795  TF1* f1ChargeGauss = new TF1("f1ChargeDist",
1796  "[0]*TMath::Gaus(x,[1],[2])",
1797  poipos_local[1] - 2 * iClusterSize,
1798  poipos_local[1] + 2 * iClusterSize);
1799  // TofChargeDistributions * fptr = new TofChargeDistributions();
1800  // TF1 * f1ChargeGauss = new TF1( "f1ChargeDist", fptr,
1801  // &TofChargeDistributions::Gauss1D,
1802  // poipos_local[1] - 2*iClusterSize,
1803  // poipos_local[1] + 2*iClusterSize,
1804  // 3
1805  // );
1806  // FIXME: wrong normalization
1807  // why is there a coefficient of 1/2 for the standard deviation?
1808  // why is the mean given in [cm] and the standard deviation dimensionless?
1809  f1ChargeGauss->SetParameters(
1810  dClustCharge / (TMath::Sqrt(2.0 * TMath::Pi() * iClusterSize / 6.0)),
1811  poipos_local[1],
1812  0.5 * iClusterSize / 3.0);
1813 
1814  double* x = new double[kiNbIntPts];
1815  double* w = new double[kiNbIntPts];
1816  f1ChargeGauss->CalcGaussLegendreSamplingPoints(kiNbIntPts, x, w, 1e-12);
1817  // Loop over strips, get the share of the charge each get, compute
1818  // the times it measure and create the digis
1819  for (Int_t iStripInd = iChannel - iNbStripsSideA;
1820  iStripInd <= iChannel + iNbStripsSideB;
1821  iStripInd++)
1822  if (0 <= iStripInd && iStripInd < iNbCh) {
1823  Int_t iCh1 = iStripInd;
1824  if (fGeoHandler->GetGeoVersion() < k14a)
1825  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
1826  CbmTofDetectorInfo xDetInfo(
1827  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
1828  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
1829  fChannelInfo = fDigiPar->GetCell(iSideChId);
1830 
1831  Double_t dStripCharge = f1ChargeGauss->IntegralFast(
1832  kiNbIntPts,
1833  x,
1834  w,
1835  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
1836  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
1837 
1838  // Fee Threshold on charge, each side gets half, each side has a gain
1839  if (fDigiBdfPar->GetFeeThreshold()
1840  < dStripCharge
1841  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc]
1842  [2 * iStripInd + 1]
1843  / 2.0) {
1844  Double_t dTimeA = dCentralTime;
1845  dTimeA +=
1846  gRandom->Gaus(0.0, fdTimeResElec)
1847 #ifdef FULL_PROPAGATION_TIME
1848  + TMath::Abs(poipos_local[0] - (+fChannelInfo->GetSizex() / 2.0))
1849 #else
1850  + (poipos_local[0])
1851 #endif
1852  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
1853  CbmTofDigi* tofDigi = new CbmTofDigi(
1854  iSM,
1855  iRpc,
1856  iChannel,
1857  dTimeA,
1858  dStripCharge
1859  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
1860  / 2.0,
1861  1,
1862  iSmType);
1863  // tofDigi->GetMatch()->AddLink(1., iPntInd);
1864  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
1865  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].push_back(
1866  data);
1867  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
1868  .push_back(iPntInd);
1869 
1870  LOG(debug) << Form(
1871  "Digimatch (%d,%d,%d,%d): size %zu, valA %d, MCt %d",
1872  iSmType,
1873  iSM,
1874  iRpc,
1875  iChannel,
1876  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
1877  .size(),
1878  iPntInd,
1879  iTrackID);
1880  } // if Side A above threshold
1881  if (fDigiBdfPar->GetFeeThreshold()
1882  < dStripCharge
1883  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iStripInd]
1884  / 2.0) {
1885  Double_t dTimeB = dCentralTime;
1886  dTimeB +=
1887  gRandom->Gaus(0.0, fdTimeResElec)
1888 #ifdef FULL_PROPAGATION_TIME
1889  + TMath::Abs(poipos_local[0] - (-fChannelInfo->GetSizex() / 2.0))
1890 #else
1891  - (poipos_local[0])
1892 #endif
1893  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
1894  CbmTofDigi* tofDigi = new CbmTofDigi(
1895  iSM,
1896  iRpc,
1897  iChannel,
1898  dTimeB,
1899  dStripCharge
1900  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel]
1901  / 2.0,
1902  0,
1903  iSmType);
1904  // tofDigi->GetMatch()->AddLink(1., iPntInd);
1905  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
1906  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iChannel].push_back(
1907  data);
1908  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel]
1909  .push_back(iPntInd);
1910  LOG(debug1) << Form(
1911  "Digimatch (%d,%d,%d,%d): size %zu, valB %d, MCt %d",
1912  iSmType,
1913  iSM,
1914  iRpc,
1915  iChannel,
1916  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
1917  .size(),
1918  iPntInd,
1919  iTrackID);
1920 
1921  } // if Side B above threshold
1922  } // if( 0 <= iStripInd && iStripInd < iNbCh )
1923  delete f1ChargeGauss;
1924  // delete []x;
1925  // delete []w;
1926  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
1927  else {
1928  // Vertical channel: A = Top and B = Bottom of the channel
1929 
1930  // Assume the width (sigma) of the gaussian is 0.5*dClusterSize/3
1931  // => 99% of the charge is in "dClusterSize"
1932  TF1* f1ChargeGauss = new TF1("f1ChargeDist",
1933  "[0]*TMath::Gaus(x,[1],[2])",
1934  poipos_local[0] - 2 * iClusterSize,
1935  poipos_local[0] + 2 * iClusterSize);
1936 
1937  // TofChargeDistributions * fptr = new TofChargeDistributions();
1938  // TF1 * f1ChargeGauss = new TF1( "f1ChargeDist", fptr,
1939  // &TofChargeDistributions::Gauss1D,
1940  // poipos_local[1] - 2*iClusterSize,
1941  // poipos_local[1] + 2*iClusterSize,
1942  // 3
1943  // );
1944  f1ChargeGauss->SetParameters(
1945  dClustCharge / (TMath::Sqrt(2.0 * TMath::Pi()) * iClusterSize / 6.0),
1946  poipos_local[0],
1947  0.5 * iClusterSize / 3.0);
1948 
1949  double* x = new double[kiNbIntPts];
1950  double* w = new double[kiNbIntPts];
1951  f1ChargeGauss->CalcGaussLegendreSamplingPoints(kiNbIntPts, x, w, 1e-12);
1952  // Loop over strips, get the share of the charge each get, compute
1953  // the times it measure and create the digis
1954  for (Int_t iStripInd = iChannel - iNbStripsSideA;
1955  iStripInd <= iChannel + iNbStripsSideB;
1956  iStripInd++)
1957  if (0 <= iStripInd && iStripInd < iNbCh) {
1958  Int_t iCh1 = iStripInd;
1959  if (fGeoHandler->GetGeoVersion() < k14a)
1960  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
1961  CbmTofDetectorInfo xDetInfo(
1962  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
1963  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
1964  fChannelInfo = fDigiPar->GetCell(iSideChId);
1965 
1966  Double_t dStripCharge = f1ChargeGauss->IntegralFast(
1967  kiNbIntPts,
1968  x,
1969  w,
1970  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
1971  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0);
1972 
1973  // Fee Threshold on charge, each side gets half, each side has a gain
1974  if (fDigiBdfPar->GetFeeThreshold()
1975  < dStripCharge
1976  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc]
1977  [2 * iStripInd + 1]
1978  / 2.0) {
1979  Double_t dTimeA = dCentralTime;
1980  dTimeA +=
1981  gRandom->Gaus(0.0, fdTimeResElec)
1982 #ifdef FULL_PROPAGATION_TIME
1983  + TMath::Abs(poipos_local[1] - (+fChannelInfo->GetSizey() / 2.0))
1984 #else
1985  + (poipos_local[1])
1986 #endif
1987  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
1988  LOG(debug) << "Create DigiA TSRC " << iSmType << iSM << iRpc
1989  << iStripInd
1990  << Form("at %f, ypos %f", dTimeA, poipos_local[1]);
1991  CbmTofDigi* tofDigi = new CbmTofDigi(
1992  iSM,
1993  iRpc,
1994  iStripInd,
1995  dTimeA,
1996  dStripCharge
1997  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iStripInd + 1]
1998  / 2.0,
1999  1,
2000  iSmType);
2001  // tofDigi->GetMatch()->AddLink(1., iPntInd);
2002  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
2003  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iStripInd + 1]
2004  .push_back(data);
2005  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iStripInd + 1]
2006  .push_back(iPntInd);
2007  LOG(debug1) << Form(
2008  "Digimatch (%d,%d,%d,%d): size %zu, valC %d, MCt %d",
2009  iSmType,
2010  iSM,
2011  iRpc,
2012  iChannel,
2013  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
2014  .size(),
2015  iPntInd,
2016  iTrackID);
2017  } // if Side A above threshold
2018  if (fDigiBdfPar->GetFeeThreshold()
2019  < dStripCharge
2020  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iStripInd]
2021  / 2.0) {
2022  Double_t dTimeB = dCentralTime;
2023  dTimeB +=
2024  gRandom->Gaus(0.0, fdTimeResElec)
2025 #ifdef FULL_PROPAGATION_TIME
2026  + TMath::Abs(poipos_local[1] - (-fChannelInfo->GetSizey() / 2.0))
2027 #else
2028  - (poipos_local[1])
2029 #endif
2030  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2031  LOG(debug) << "Create DigiB TSRC " << iSmType << iSM << iRpc
2032  << iStripInd
2033  << Form("at %f, ypos %f", dTimeB, poipos_local[1]);
2034  CbmTofDigi* tofDigi = new CbmTofDigi(
2035  iSM,
2036  iRpc,
2037  iStripInd,
2038  dTimeB,
2039  dStripCharge
2040  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iStripInd]
2041  / 2.0,
2042  0,
2043  iSmType);
2044  // tofDigi->GetMatch()->AddLink(1., iPntInd);
2045  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
2046  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iStripInd].push_back(
2047  data);
2048  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iStripInd]
2049  .push_back(iPntInd);
2050  LOG(debug1) << Form(
2051  "Digimatch (%d,%d,%d,%d): size %zu, valE %d, MCt %d",
2052  iSmType,
2053  iSM,
2054  iRpc,
2055  iChannel,
2056  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
2057  .size(),
2058  iPntInd,
2059  iTrackID);
2060  } // if Side B above threshold
2061  } // if( 0 <= iStripInd && iStripInd < iNbCh )
2062  delete f1ChargeGauss;
2063  delete[] x;
2064  delete[] w;
2065  } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2066  } // for (Int_t iPntInd = 0; iPntInd < nTofPoint; iPntInd++ )
2067 
2068  // Clear the Track to channel temporary storage
2069  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
2070  for (UInt_t uTrkInd = 0; uTrkInd < fvlTrckChAddr.size(); uTrkInd++) {
2071  fvlTrckChAddr[uTrkInd].clear();
2072  fvlTrckRpcAddr[uTrkInd].clear();
2073  fvlTrckRpcTime[uTrkInd].clear();
2074  }
2075  fvlTrckChAddr.clear();
2076  fvlTrckRpcAddr.clear();
2077  fvlTrckRpcTime.clear();
2078  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
2079 
2080  return kTRUE;
2081 }
2082 /************************************************************************************/
2083 // Functions for a simple "Flat disc" cluster approximation
2084 Bool_t CbmTofDigitize::DigitizeFlatDisc() {
2085  // Uniform distribution in ]0;x]
2086  // gRandom->Uniform(x);
2087  // Gauss distribution in of mean m, sigma s
2088  // gRandom->Gauss(m, s);
2089 
2090  CbmTofPoint* pPoint;
2091  CbmMCTrack* pMcTrk;
2092 
2093  Int_t nTofPoint = fTofPointsColl->GetEntries();
2094  Int_t nMcTracks = fMcTracksColl->GetEntries();
2095 
2096  // Prepare the temporary storing of the Track/Point/Digi info
2097  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
2098  fvlTrckChAddr.resize(nMcTracks);
2099  fvlTrckRpcAddr.resize(nMcTracks);
2100  fvlTrckRpcTime.resize(nMcTracks);
2101  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
2102  for (Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++) {
2103  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
2104  fvlTrckChAddr[iTrkInd].clear();
2105  fvlTrckRpcAddr[iTrkInd].clear();
2106  fvlTrckRpcTime[iTrkInd].clear();
2107  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
2108  } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++)
2109 
2110  if (fbMcTrkMonitor) {
2111  // Debug Info printout
2112  Int_t iNbTofTracks = 0;
2113  Int_t iNbTofTracksPrim = 0;
2114 
2115  LOG(debug1) << "CbmTofDigitize::DigitizeFlatDisc: " << nTofPoint
2116  << " points in Tof for this event with " << nMcTracks
2117  << " MC tracks ";
2118  for (Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++) {
2119  pMcTrk = (CbmMCTrack*) fMcTracksColl->At(iTrkInd);
2120  if (0 < pMcTrk->GetNPoints(ECbmModuleId::kTof)) iNbTofTracks++;
2121  if (0 < pMcTrk->GetNPoints(ECbmModuleId::kTof)
2122  && -1 == pMcTrk->GetMotherId())
2123  iNbTofTracksPrim++;
2124  } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++)
2125 
2126  //Some numbers on TOF distributions
2127  LOG(debug1) << "CbmTofDigitize::DigitizeFlatDisc: " << iNbTofTracks
2128  << " tracks in Tof ";
2129  LOG(debug1) << "CbmTofDigitize::DigitizeFlatDisc: " << iNbTofTracksPrim
2130  << " tracks in Tof from vertex";
2131  } // if( fbMcTrkMonitor )
2132 
2133  // Tof Point Info
2134  Int_t iDetId, iSmType, iSM, iRpc, iChannel, iGap;
2135  Int_t iTrackID, iChanId;
2136  TVector3 vPntPos;
2137  // Cluster Info
2138  Double_t dClusterSize;
2139  Double_t dClustCharge;
2140  // Digi
2141  // CbmTofDigi * pTofDigi; // <- Compressed digi (64 bits)
2142  // CbmTofDigi * pTofDigiExp; // <- Expanded digi
2143  // General geometry info
2144  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
2145  Int_t iNbSm, iNbRpc, iNbCh;
2146  Int_t iChType;
2147 
2148  // Start jitter: Same contrib. for all points in same event
2149  // FIXME
2150  // Actually, event start time reconstruction should not be done in the ToF
2151  // digitizer, neither in event-based nor in time-based mode. The timestamp
2152  // of a CbmTofDigi object is not the time difference between signals in a
2153  // (hypothetical) start detector system and the MRPC wall but rather
2154  // represents the detection point in time in an MRPC measured with respect to
2155  // the start of the run.
2156  // For compatibility reasons, we continue adding a start detector jitter to
2157  // digi timestamps in event-based mode. In time-based mode, event start time
2158  // reconstruction is not considered here.
2159  Double_t dStartJitter = 0.;
2160  if (!fbTimeBasedOutput) {
2161  dStartJitter = gRandom->Gaus(0.0, fDigiBdfPar->GetStartTimeRes());
2162  }
2163 
2164  for (Int_t iPntInd = 0; iPntInd < nTofPoint; iPntInd++) {
2165  // Get a pointer to the TOF point
2166  pPoint = (CbmTofPoint*) fTofPointsColl->At(iPntInd);
2167  if (NULL == pPoint) {
2168  LOG(WARNING) << "CbmTofDigitize::DigitizeFlatDisc => Be careful: hole in "
2169  "the CbmTofPoint TClonesArray!"
2170  << endl;
2171  continue;
2172  } // if( pPoint )
2173 
2174  // Get all channel info
2175  iDetId = pPoint->GetDetectorID();
2176  iSmType = fGeoHandler->GetSMType(iDetId);
2177  iRpc = fGeoHandler->GetCounter(iDetId);
2178  iChannel = fGeoHandler->GetCell(iDetId);
2179  if (fGeoHandler->GetGeoVersion() < k14a)
2180  iChannel--; // Again, channel going from 1 to nbCh instead of 0 to nbCh - 1 ?!?!?
2181  iGap = fGeoHandler->GetGap(iDetId);
2182  iSM = fGeoHandler->GetSModule(iDetId);
2183  iChanId = fGeoHandler->GetCellId(iDetId);
2184  iTrackID = pPoint->GetTrackID();
2185 
2186  iNbSm = fDigiBdfPar->GetNbSm(iSmType);
2187  iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
2188  iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
2189  iChType = fDigiBdfPar->GetChanType(iSmType, iRpc);
2190  Int_t iTrkId = pPoint->GetTrackID();
2191 
2192  // Catch case where the MC track was removed from the transport stack
2193  if (iTrkId < 0) {
2194  if (fbAllowPointsWithoutTrack)
2195  continue;
2196  else
2197  LOG(fatal)
2198  << "CbmTofDigitize::DigitizeDirectClusterSize => TofPoint without "
2199  "valid MC track Index, "
2200  << " track was probably cut at the transport level! "
2201  << " Pnt Idx: " << iPntInd << " Trk Idx: " << iTrkId << "\n"
2202  << "=============> To allow this kind of points and simply jump "
2203  "them, "
2204  << "call CbmTofDigitize::AllowPointsWithoutTrack() in your macro!!";
2205  } // if( iTrkId < 0 )
2206 
2207  if (fbMcTrkMonitor) {
2208  // Get pointer to the MC-Track info
2209  pMcTrk = (CbmMCTrack*) fMcTracksColl->At(iTrkId);
2210 
2211  LOG(debug1) << Form("CbmTofDigitize => det ID 0x%08x", iDetId)
2212  << ", GeoVersion " << fGeoHandler->GetGeoVersion()
2213  << ", SMType: " << iSmType << " SModule: " << iSM << " of "
2214  << iNbSm + 1 << " Module: " << iRpc << " of " << iNbRpc + 1
2215  << " Gap: " << iGap << " Cell: " << iChannel << " of "
2216  << iNbCh + 1;
2217 
2218  // Jump all tracks not making 8 points for test
2219  // if( 8 != pMcTrk->GetNPoints(ECbmModuleId::kTof) )
2220  // continue;
2221 
2222  // Keep only primaries
2223  if (kTRUE == fDigiBdfPar->UseOnlyPrimaries())
2224  if (-1 != pMcTrk->GetMotherId()) continue;
2225  } // if( fbMcTrkMonitor )
2226 
2227  if (iSmType < 0. || iNbSmTypes <= iSmType || iSM < 0. || iNbSm <= iSM
2228  || iRpc < 0. || iNbRpc <= iRpc || iChannel < 0. || iNbCh <= iChannel) {
2229  LOG(error) << Form("CbmTofDigitize => det ID 0x%08x", iDetId)
2230  << ", GeoVersion " << fGeoHandler->GetGeoVersion()
2231  << ", SMType: " << iSmType;
2232  LOG(error) << " SModule: " << iSM << " of " << iNbSm + 1;
2233  LOG(error) << " Module: " << iRpc << " of " << iNbRpc + 1;
2234  LOG(error) << " Gap: " << iGap;
2235  LOG(fatal) << " Cell: " << iChannel << " of " << iNbCh + 1;
2236  continue;
2237  } // if error on channel data
2238  LOG(debug2) << "CbmTofDigitize::DigitizeFlatDisc => Check Point " << iPntInd
2239  << " Sm type " << iSmType << " SM " << iSM << " Rpc " << iRpc
2240  << " Channel " << iChannel << " Type " << iChType << " Orient. "
2241  << fDigiBdfPar->GetChanOrient(iSmType, iRpc);
2242 
2243  // Check if there was already a Digi from the same track created
2244  // with this channel as main channel
2245  ULong64_t uAddr =
2246  CbmTofAddress::GetUniqueAddress(iSM, iRpc, iChannel, 0, iSmType);
2247  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
2248  Bool_t bFoundIt = kFALSE;
2249  // Check if this track ID is bigger than size of McTracks TClonesArray (maybe happen with PSD cut!!)
2250  // If yes, resize the array to reach the appropriate number of fields
2251  // Cast of TrackId is safe as we already check for "< 0" case
2252  // TODO: Is this check still required when Id < 0 are jumped?
2253  if (fvlTrckChAddr.size() <= static_cast<UInt_t>(iTrkId))
2254  fvlTrckChAddr.resize(iTrkId + 1);
2255 
2256  for (UInt_t uTrkMainCh = 0; uTrkMainCh < fvlTrckChAddr[iTrkId].size();
2257  uTrkMainCh++)
2258  if (uAddr == fvlTrckChAddr[iTrkId][uTrkMainCh]) {
2259  bFoundIt = kTRUE;
2260  break;
2261  }
2262  // If it is the case, no need to redo the full stuff for this gap
2263  if (kTRUE == bFoundIt) continue;
2264  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
2265 
2266  // Probability that the RPC detect the hit
2267  // if( fDigiBdfPar->GetEfficiency(iSmType) < gRandom->Uniform(1) )
2268  // continue;
2269 
2270  // Probability that the gap detect the hit
2271  // For now use a simple gap treatment (cf CbmTofDigiBdfPar):
2272  // - a single fired gap fires the channel
2273  // - all gaps have exactly equivalent efficiency/firing probability
2274  // - independent gap firing (no charge accumulation or such)
2275  // The probability for a hit not to fire at all is then
2276  // (1-p)^NGaps with p = gap efficiency and Ngaps the number of gaps in this RPC
2277  // This results in the relation: p = 1 - (1 - P)^(1/Ngaps)
2278  // with P = RPC efficiency from beamtime
2279  LOG(debug2) << "CbmTofDigitize::DigitizeFlatDisc => Eff = "
2280  << fDigiBdfPar->GetGapEfficiency(iSmType, iRpc);
2281 
2282  if (fDigiBdfPar->GetGapEfficiency(iSmType, iRpc) < gRandom->Uniform(1))
2283  continue;
2284 
2285  // Save the address in vector so that this channel is not redone later for the
2286  // eventual other gaps TofPoint
2287  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk())
2288  fvlTrckChAddr[iTrkId].push_back(uAddr);
2289 
2290  // Get TofPoint Position
2291  pPoint->Position(vPntPos);
2292  fChannelInfo = fDigiPar->GetCell(iChanId);
2293  if (NULL == fChannelInfo) {
2294  LOG(WARNING)
2295  << "CbmTofDigitize::DigitizeFlatDisc: No DigPar for iChanId = "
2296  << Form("0x%08x, addr 0x%08x", iChanId, (unsigned int) iDetId);
2297  continue;
2298  }
2299  // Master -> Local
2300  TGeoNode* fNode = // prepare global->local trafo
2301  gGeoManager->FindNode(
2302  fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
2303  LOG(debug1) << Form(" TofDigitizerBDF:: (%3d,%3d,%3d,%3d) - node at "
2304  "(%6.1f,%6.1f,%6.1f) : %p Pos(%6.1f,%6.1f,%6.1f)",
2305  iSmType,
2306  iSM,
2307  iRpc,
2308  iChannel,
2309  fChannelInfo->GetX(),
2310  fChannelInfo->GetY(),
2311  fChannelInfo->GetZ(),
2312  fNode,
2313  vPntPos.X(),
2314  vPntPos.Y(),
2315  vPntPos.Z());
2316  gGeoManager->GetCurrentNode();
2317  gGeoManager->GetCurrentMatrix();
2318 
2319  Double_t refpos[3] = {
2320  fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ()};
2321  Double_t refpos_local[3];
2322  gGeoManager->MasterToLocal(refpos, refpos_local);
2323  Double_t poipos[3] = {vPntPos.X(), vPntPos.Y(), vPntPos.Z()};
2324  Double_t poipos_local[3];
2325  gGeoManager->MasterToLocal(poipos, poipos_local);
2326  LOG(debug1) << Form(" TofDigitizerBDF:: DetId 0x%08x, ChanId 0x%08x, local "
2327  "pos (%5.1f,%5.1f,%5.1f) refpos (%5.1f,%5.1f,%5.1f) ",
2328  iDetId,
2329  iChanId,
2330  poipos_local[0],
2331  poipos_local[1],
2332  poipos_local[2],
2333  refpos_local[0],
2334  refpos_local[1],
2335  refpos_local[2]);
2336  for (Int_t i = 0; i < 3; i++)
2337  poipos_local[i] -=
2338  refpos_local[i]; // correct for transformation failure, FIXME!
2339 
2340  // Generate a Cluster radius in [cm]
2341  dClusterSize = GenerateClusterRadius(iSmType, iRpc);
2342  // Cut on the higher radius because otherwise the big clusters at the edge
2343  // of the counter generate hits totally off (because the main strip do not get
2344  // more charge than the side strips anymore). Cut value fully empirical.
2345  while (dClusterSize < 0.0001 || 3.0 < dClusterSize)
2346  // Should not happen without error message
2347  // But Landau can give really small values
2348  dClusterSize = GenerateClusterRadius(iSmType, iRpc);
2349 
2350  // Generate Cluster charge as ToT[ps]
2351  dClustCharge = fh1ClusterTotProb[iSmType]->GetBinContent(
2352  fh1ClusterTotProb[iSmType]->FindBin(gRandom->Uniform(1)));
2353 
2354  //Total cluster area (used to calculate the charge fraction in each channel)
2355  Double_t dClustArea = dClusterSize * dClusterSize * TMath::Pi();
2356 
2357  // Calculate the fraction of the charge in central channel
2358  Double_t dChargeCentral =
2359  dClustCharge
2360  * ComputeClusterAreaOnChannel(
2361  iChanId, dClusterSize, poipos_local[0], poipos_local[1]);
2362  LOG(debug2) << "CbmTofDigitize::DigitizeFlatDisc: ChargeCentral "
2363  << dChargeCentral << ", " << dClustCharge
2364  << Form(", 0x%08x", iChanId) << ", " << dClusterSize << ", "
2365  << poipos_local[0] << ", " << poipos_local[1] << ", "
2366  << poipos_local[2];
2367 
2368  dChargeCentral /= dClustArea;
2369  if (dClustCharge + 0.0000001 < dChargeCentral) {
2370  LOG(error) << "CbmTofDigitize::DigitizeFlatDisc => Central Charge "
2371  << dChargeCentral << " " << dClustCharge << " "
2372  << dClustCharge - dChargeCentral << " "
2373  << (ComputeClusterAreaOnChannel(
2374  iChanId, dClusterSize, poipos_local[0], poipos_local[1]))
2375  << " " << dClustArea;
2376  LOG(error) << "CbmTofDigitize::DigitizeFlatDisc => Check Point "
2377  << iPntInd << " Sm type " << iSmType << " SM " << iSM
2378  << " Rpc " << iRpc << " Channel " << iChannel << " Type "
2379  << iChType << " Orient. "
2380  << fDigiBdfPar->GetChanOrient(iSmType, iRpc);
2381  } // if( dClustCharge +0.0000001 < dChargeCentral )
2382 
2383  // Calculate the time for the central channel
2384  Double_t dCentralTime;
2385  // Check if there was already a Digi from the same track created in this RPC
2386  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
2387  ULong64_t uRpcAddr =
2388  CbmTofAddress::GetUniqueAddress(iSM, iRpc, 0, 0, iSmType);
2389  Bool_t bFoundIt = kFALSE;
2390  UInt_t uTrkRpcPair = 0;
2391  for (uTrkRpcPair = 0; uTrkRpcPair < fvlTrckRpcAddr[iTrkId].size();
2392  uTrkRpcPair++)
2393  if (uRpcAddr == fvlTrckRpcAddr[iTrkId][uTrkRpcPair]) {
2394  bFoundIt = kTRUE;
2395  break;
2396  }
2397  // If it is the case, we should reuse the timing already assigned to this track
2398  if (kTRUE == bFoundIt)
2399  dCentralTime = fvlTrckRpcTime[iTrkId][uTrkRpcPair];
2400  else {
2401  dCentralTime =
2402  pPoint->GetTime()
2403  + gRandom->Gaus(0.0, fDigiBdfPar->GetResolution(iSmType))
2404  + dStartJitter; // Same contrib. for all points in same event
2405  fvlTrckRpcAddr[iTrkId].push_back(uRpcAddr);
2406  fvlTrckRpcTime[iTrkId].push_back(dCentralTime);
2407  } // No Digi yet in this RPC for this Track
2408  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
2409  else
2410  dCentralTime =
2411  pPoint->GetTime()
2412  + gRandom->Gaus(0.0, fDigiBdfPar->GetResolution(iSmType))
2413  + dStartJitter; // Same contrib. for all points in same event
2414 
2415  // FIXME: not sure if this limit does not destroy rate estimates and such
2416  // if( 1e6 < dCentralTime )
2417  // continue;
2418 
2419  // Calculate propagation time(s) to the readout point(s)
2420  if (0 == iChType) {
2421  // Strips (readout on both side)
2422  // Assume that the bottom/left strip have lower channel index
2423  // E.g: ... | i | i+1 | ...
2424  // or ... y
2425  // ----- ^
2426  // i+1 |
2427  // ----- --> x
2428  // i
2429  // -----
2430  // ...
2431 
2432  Double_t dTimeA = dCentralTime;
2433  Double_t dTimeB = dCentralTime;
2434  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
2435  // Horizontal channel: A = Right and B = Left of the channel
2436  dTimeA +=
2437  gRandom->Gaus(0.0, fdTimeResElec)
2438 #ifdef FULL_PROPAGATION_TIME
2439  + TMath::Abs(poipos_local[0] - (+fChannelInfo->GetSizex() / 2.0))
2440 #else
2441  + (poipos_local[0])
2442 #endif
2443  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2444  dTimeB +=
2445  gRandom->Gaus(0.0, fdTimeResElec)
2446 #ifdef FULL_PROPAGATION_TIME
2447  + TMath::Abs(poipos_local[0] - (-fChannelInfo->GetSizex() / 2.0))
2448 #else
2449  - (poipos_local[0])
2450 #endif
2451  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2452  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2453  else {
2454  // Vertical channel: A = Top and B = Bottom of the channel
2455  dTimeA +=
2456  gRandom->Gaus(0.0, fdTimeResElec)
2457 #ifdef FULL_PROPAGATION_TIME
2458  + TMath::Abs(poipos_local[1] - (+fChannelInfo->GetSizey() / 2.0))
2459 #else
2460  // + ( poipos_local[1] - fChannelInfo->GetY() )
2461  + poipos_local[1]
2462 #endif
2463  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2464  dTimeB +=
2465  gRandom->Gaus(0.0, fdTimeResElec)
2466 #ifdef FULL_PROPAGATION_TIME
2467  + TMath::Abs(poipos_local[1] - (-fChannelInfo->GetSizey() / 2.0))
2468 #else
2469  // - ( poipos_local[1] - fChannelInfo->GetY() )
2470  - poipos_local[1]
2471 #endif
2472  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2473  } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2474  if (fDigiBdfPar->GetFeeThreshold()
2475  <= dChargeCentral
2476  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
2477  / 2.0) {
2478  CbmTofDigi* tofDigi = new CbmTofDigi(
2479  iSM,
2480  iRpc,
2481  iChannel,
2482  dTimeA,
2483  dChargeCentral
2484  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
2485  / 2.0,
2486  1,
2487  iSmType);
2488  //nh,crashes
2489  // tofDigi->GetMatch()->AddLink(1., iPntInd);
2490  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
2491  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].push_back(
2492  data);
2493  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
2494  .push_back(iPntInd);
2495  LOG(debug1) << Form(
2496  "Digimatch (%d,%d,%d,%d): size %zu, MCp %d, MCt %d, TA %6.2f, TotA "
2497  "%6.1f, Ypos %6.1f, Vsig %6.1f ",
2498  iSmType,
2499  iSM,
2500  iRpc,
2501  iChannel,
2502  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].size(),
2503  iPntInd,
2504  iTrackID,
2505  dTimeA,
2506  dChargeCentral
2507  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
2508  / 2.0,
2509  poipos_local[1],
2510  fvdSignalVelocityRpc[iSmType][iSM][iRpc]);
2511 
2512  } // charge ok, TimeA
2513  if (fDigiBdfPar->GetFeeThreshold()
2514  <= dChargeCentral
2515  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel]
2516  / 2.0) {
2517  CbmTofDigi* tofDigi = new CbmTofDigi(
2518  iSM,
2519  iRpc,
2520  iChannel,
2521  dTimeB,
2522  dChargeCentral
2523  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel] / 2.0,
2524  0,
2525  iSmType);
2526  // tofDigi->GetMatch()->AddLink(1., iPntInd);
2527  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
2528  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iChannel].push_back(data);
2529  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel].push_back(
2530  iPntInd);
2531  LOG(debug1) << Form(
2532  "Digimatch (%d,%d,%d,%d): size %zu, MCp %d, MCt %d, TB %6.2f, TotB "
2533  "%6.1f",
2534  iSmType,
2535  iSM,
2536  iRpc,
2537  iChannel,
2538  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].size(),
2539  iPntInd,
2540  iTrackID,
2541  dTimeB,
2542  dChargeCentral
2543  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
2544  / 2.0);
2545  } // charge ok, TimeB
2546  } // if( 0 == iChType)
2547  else {
2548  // Pad (Single side readout)
2549  // Assume that the bottom/left pads have lower channel index
2550  // E.g: for a 2*4 pads counter, pads 0-3 are the bottom/left ones and 4-7 the
2551  // top/right one with reversed numbering:
2552  // -----------------
2553  // row 1 | 7 | 6 | 5 | 4 | y
2554  // ----------------- ^
2555  // row 0 | 0 | 1 | 2 | 3 | |
2556  // ----------------- --> x
2557  // or ---------
2558  // | 4 | 3 |
2559  // ---------
2560  // | 5 | 2 |
2561  // ---------
2562  // | 6 | 1 |
2563  // ---------
2564  // | 7 | 0 |
2565  // ---------
2566  // row 1 0
2567  // Also assume that the readout happens at the middle of the readout edge
2568 
2569  Double_t dClustToReadout = 0.0;
2570  Double_t dPadTime = dCentralTime;
2571  // First calculate the propagation time to the center of the pad base
2572  if (iChannel < iNbCh / 2.0) {
2573  // First row
2574  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
2575  // Vertical => base = right edge
2576  dClustToReadout = TMath::Sqrt(
2577  TMath::Power(poipos_local[1], 2)
2578  + TMath::Power(poipos_local[0] - (+fChannelInfo->GetSizex() / 2.0),
2579  2));
2580  else // Horizontal => base = bottom edge
2581  dClustToReadout = TMath::Sqrt(
2582  TMath::Power(poipos_local[0], 2)
2583  + TMath::Power(poipos_local[1] - (-fChannelInfo->GetSizey() / 2.0),
2584  2));
2585  } // if( iChannel < iNbCh/2.0 )
2586  else {
2587  // Second row
2588  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
2589  // Vertical => base = left edge
2590  dClustToReadout = TMath::Sqrt(
2591  TMath::Power(poipos_local[1], 2)
2592  + TMath::Power(poipos_local[0] - (-fChannelInfo->GetSizex() / 2.0),
2593  2));
2594  else // Horizontal => base = upper edge
2595  dClustToReadout = TMath::Sqrt(
2596  TMath::Power(poipos_local[0], 2)
2597  + TMath::Power(poipos_local[1] - (+fChannelInfo->GetSizey() / 2.0),
2598  2));
2599  } // else of if( iChannel < iNbCh/2.0 )
2600 
2601  dPadTime += gRandom->Gaus(0.0, fdTimeResElec)
2602  + dClustToReadout / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2603 
2604  if (fDigiBdfPar->GetFeeThreshold()
2605  <= dChargeCentral
2606  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iChannel]) {
2607  CbmTofDigi* tofDigi = new CbmTofDigi(
2608  iSM,
2609  iRpc,
2610  iChannel,
2611  dPadTime,
2612  dChargeCentral
2613  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iChannel],
2614  0,
2615  iSmType);
2616  // tofDigi->GetMatch()->AddLink(1., iPntInd);
2617  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
2618  fStorDigi[iSmType][iSM * iNbRpc + iRpc][iChannel].push_back(data);
2619  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iChannel].push_back(
2620  iPntInd);
2621  LOG(debug1) << Form(
2622  "Digimatch (%d,%d,%d,%d): MCpi %zu, valE %d, MCti %d",
2623  iSmType,
2624  iSM,
2625  iRpc,
2626  iChannel,
2627  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].size(),
2628  iPntInd,
2629  iTrackID);
2630 
2631  } // charge ok
2632  } // else of if( 0 == iChType)
2633 
2634  // Loop over neighboring channel to find if they have overlap with
2635  // the cluster ( and if yes which fraction of the total charge they
2636  // got)
2637  if (0 == iChType) {
2638  // Strips (readout on both side)
2639  // Loop in decreasing channel index until a channel with right/upper edge farther
2640  // from cluster center than cluster radius is found
2641  // Loop in increasing channel index until a channel with left/lower edge farther
2642  // from cluster center than cluster radius is found
2643 
2644  Int_t iMinChanInd = iChannel - 1;
2645  Int_t iMaxChanInd = iChannel + 1;
2646 
2647  Double_t dClusterDist = 0.0;
2648  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
2649  // Horizontal channels: First go down, then up
2650  while (0 <= iMinChanInd) {
2651  dClusterDist = TMath::Abs(
2652  poipos_local[1]
2653  - (+fChannelInfo->GetSizey() * (iMinChanInd - iChannel + 0.5)));
2654  if (dClusterDist < dClusterSize)
2655  iMinChanInd--;
2656  else
2657  break;
2658  }
2659  while (iMaxChanInd < iNbCh) {
2660  dClusterDist = TMath::Abs(
2661  poipos_local[1]
2662  - (+fChannelInfo->GetSizey() * (iMaxChanInd - iChannel - 0.5)));
2663  if (dClusterDist < dClusterSize)
2664  iMaxChanInd++;
2665  else
2666  break;
2667  }
2668  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2669  else {
2670  // Vertical channels: First go to the left, then to the right
2671  while (0 <= iMinChanInd) {
2672  dClusterDist = TMath::Abs(
2673  poipos_local[0]
2674  - (+fChannelInfo->GetSizex() * (iMinChanInd - iChannel + 0.5)));
2675  if (dClusterDist < dClusterSize)
2676  iMinChanInd--;
2677  else
2678  break;
2679  }
2680  while (iMaxChanInd < iNbCh) {
2681  dClusterDist = TMath::Abs(
2682  poipos_local[0]
2683  - (+fChannelInfo->GetSizex() * (iMaxChanInd - iChannel - 0.5)));
2684  if (dClusterDist < dClusterSize)
2685  iMaxChanInd++;
2686  else
2687  break;
2688  }
2689  } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2690 
2691  // Loop over all channels inside the interval ]iMinChanInd;iMaxChanInd[ except central channel,
2692  // using the overlap area with cluster to get the charge and adding a Digi for all channels
2693  for (Int_t iChanInd = iMinChanInd + 1; iChanInd < iMaxChanInd;
2694  iChanInd++) {
2695  if (iChanInd == iChannel) continue;
2696 
2697  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
2698  // ... don't ask me why ... Reason found in TofMC and TofGeoHandler
2699  // CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iChanInd + 1);
2700  Int_t iCh1 = iChanInd;
2701  if (fGeoHandler->GetGeoVersion() < k14a)
2702  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
2703  CbmTofDetectorInfo xDetInfo(
2704  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
2705  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
2706 
2707  fChannelInfo = fDigiPar->GetCell(iSideChId);
2708 
2709  if (NULL == fChannelInfo) {
2710  LOG(error) << "CbmTofDigitize::DigitizeFlatDisc: Invalid channel "
2711  << iSideChId << " " << iSmType << " " << iSM << " " << iRpc
2712  << " " << iChanInd + 1;
2713  continue;
2714  }
2715 
2716  // Calculate the fraction of the charge in this channel
2717  Double_t dChargeSideCh =
2718  dClustCharge
2719  * ComputeClusterAreaOnChannel(
2720  iSideChId, dClusterSize, poipos_local[0], poipos_local[1]);
2721  dChargeSideCh /= dClustArea;
2722  if (dClustCharge + 0.0000001 < dChargeSideCh) {
2723  LOG(error) << "CbmTofDigitize::DigitizeFlatDisc => Side Charge "
2724  << dChargeSideCh << " " << dClustCharge << " "
2725  << dClustArea;
2726  LOG(error) << "CbmTofDigitize::DigitizeFlatDisc => Check Point "
2727  << iPntInd << " Sm type " << iSmType << " SM " << iSM
2728  << " Rpc " << iRpc << " Channel " << iChanInd << " Type "
2729  << iChType << " Orient. "
2730  << fDigiBdfPar->GetChanOrient(iSmType, iRpc);
2731  } // if( dClustCharge + 0.0000001 < dChargeSideCh )
2732 
2733  // Strips = readout on both side
2734  Double_t dTimeA = dCentralTime;
2735  Double_t dTimeB = dCentralTime;
2736 
2737  LOG(debug1) << Form(" TofDigitizerBDF:: chrg %7.1f, times %5.1f,%5.1f, "
2738  "pos %5.1f,%5.1f,%5.1f ",
2739  dChargeCentral,
2740  dTimeA,
2741  dTimeB,
2742  poipos_local[0],
2743  poipos_local[1],
2744  poipos_local[2]);
2745 
2746  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
2747  // Horizontal channel: A = Right and B = Left of the channel
2748  dTimeA +=
2749  gRandom->Gaus(0.0, fdTimeResElec)
2750 #ifdef FULL_PROPAGATION_TIME
2751  + TMath::Abs(poipos_local[0] - (+fChannelInfo->GetSizex() / 2.0))
2752 #else
2753  + (poipos_local[0])
2754 #endif
2755  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2756  dTimeB +=
2757  gRandom->Gaus(0.0, fdTimeResElec)
2758 #ifdef FULL_PROPAGATION_TIME
2759  + TMath::Abs(poipos_local[0] - (-fChannelInfo->GetSizex() / 2.0))
2760 #else
2761  - (poipos_local[0])
2762 #endif
2763  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2764  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2765  else {
2766 
2767  // Vertical channel: A = Top and B = Bottom of the channel
2768  dTimeA +=
2769  gRandom->Gaus(0.0, fdTimeResElec)
2770 #ifdef FULL_PROPAGATION_TIME
2771  + TMath::Abs(poipos_local[1] - (+fChannelInfo->GetSizey() / 2.0))
2772 #else
2773  // + ( poipos_local[1] - fChannelInfo->GetY() )
2774  + poipos_local[1]
2775 #endif
2776  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2777  dTimeB +=
2778  gRandom->Gaus(0.0, fdTimeResElec)
2779 #ifdef FULL_PROPAGATION_TIME
2780  + TMath::Abs(poipos_local[1] - (-fChannelInfo->GetSizey() / 2.0))
2781 #else
2782  // - ( poipos_local[1] - fChannelInfo->GetY() )
2783  - poipos_local[1]
2784 #endif
2785  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
2786  } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2787  LOG(debug1) << Form(
2788  " TofDigitizerBDF:: chrg %7.1f, gain %7.1f, thr %7.1f times "
2789  "%5.1f,%5.1f ",
2790  dChargeCentral,
2791  fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1],
2792  fDigiBdfPar->GetFeeThreshold(),
2793  dTimeA,
2794  dTimeB);
2795  // Fee Threshold on charge
2796  if (fDigiBdfPar->GetFeeThreshold()
2797  <= dChargeSideCh
2798  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd + 1]
2799  / 2.0) {
2800  CbmTofDigi* tofDigi = new CbmTofDigi(
2801  iSM,
2802  iRpc,
2803  iChanInd,
2804  dTimeA,
2805  dChargeSideCh
2806  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd + 1]
2807  / 2.0,
2808  1,
2809  iSmType);
2810  // tofDigi->GetMatch()->AddLink(1., iPntInd);
2811  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
2812  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd + 1].push_back(
2813  data);
2814  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd + 1]
2815  .push_back(iPntInd);
2816  LOG(debug1) << Form(
2817  "Digimatch (%d,%d,%d,%d): size %zu, MCpA %d, MCtt %d",
2818  iSmType,
2819  iSM,
2820  iRpc,
2821  iChanInd,
2822  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd + 1]
2823  .size(),
2824  iPntInd,
2825  iTrackID);
2826 
2827  } // if( charge above threshold)
2828  if (fDigiBdfPar->GetFeeThreshold()
2829  <= dChargeSideCh
2830  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd]
2831  / 2.0) {
2832  CbmTofDigi* tofDigi = new CbmTofDigi(
2833  iSM,
2834  iRpc,
2835  iChanInd,
2836  dTimeB,
2837  dChargeSideCh
2838  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd] / 2.0,
2839  0,
2840  iSmType);
2841  // tofDigi->GetMatch()->AddLink(1., iPntInd);
2842  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
2843  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd].push_back(data);
2844  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd].push_back(
2845  iPntInd);
2846  LOG(debug1) << Form(
2847  "Digimatch (%d,%d,%d,%d): size %zu, MCpB %d, MCtB %d",
2848  iSmType,
2849  iSM,
2850  iRpc,
2851  iChanInd,
2852  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChanInd + 1]
2853  .size(),
2854  iPntInd,
2855  iTrackID);
2856 
2857  } // if( charge above threshold)
2858  } // for( Int_t iChanInd = iMinChanInd + 1; iChanInd < iMaxChanInd; iChanInd++ )
2859  } // if( 0 == iChType)
2860  else {
2861  // Pad (Single side readout)
2862  // Assume that the bottom/left pads have lower channel index
2863  // E.g: for a 2*4 pads counter, pads 0-3 are the bottom/left ones and 4-7 the
2864  // top/right one with reversed numbering
2865  // -----------------
2866  // row 1 | 7 | 6 | 5 | 4 | y
2867  // ----------------- ^
2868  // row 0 | 0 | 1 | 2 | 3 | |
2869  // ----------------- --> x
2870  // or ---------
2871  // | 4 | 3 |
2872  // ---------
2873  // | 5 | 2 |
2874  // ---------
2875  // | 6 | 1 |
2876  // ---------
2877  // | 7 | 0 |
2878  // ---------
2879  // row 1 0
2880 
2881  // Loop in decreasing channel index in same row, find min = 1st 0 channel
2882  // Loop in increasing channel index in same row, find max = 1st 0 channel
2883  Int_t iMinChanInd = iChannel;
2884  Int_t iMaxChanInd = iChannel;
2885 
2886  Double_t dClusterDist = 0;
2887  Int_t iRow;
2888  Bool_t bCheckOtherRow = kFALSE;
2889  if (iChannel < iNbCh / 2.0)
2890  iRow = 0;
2891  else
2892  iRow = 1;
2893 
2894  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
2895  // Vertical => base = right edge for row 0 and left edge for row 1
2896  while (dClusterDist < dClusterSize
2897  && iRow * iNbCh / 2.0 < iMinChanInd) {
2898  iMinChanInd--;
2899  dClusterDist =
2900  TMath::Abs(poipos_local[1]
2901  - (+fChannelInfo->GetSizey()
2902  * (iMinChanInd - iChannel + (1 - 2 * iRow) / 2)));
2903  } // while upper/lower edge inside cluster and index not out of rpc
2904  dClusterDist = 0;
2905  while (dClusterDist < dClusterSize
2906  && iMaxChanInd < (1 + iRow) * iNbCh / 2.0 - 1) {
2907  iMaxChanInd++;
2908  dClusterDist =
2909  TMath::Abs(poipos_local[1]
2910  - (+fChannelInfo->GetSizey()
2911  * (iMaxChanInd - iChannel - (1 - 2 * iRow) / 2)));
2912  } // while lower/upper edge inside cluster and index not out of rpc
2913 
2914  // Test if Pad in diff row but same column as central pad has cluster overlap
2915  // if Yes => Loop from min+1 to max-1 equivalents in the opposite row
2916  dClusterDist = TMath::Abs(
2917  poipos_local[0] - (+fChannelInfo->GetSizex() * (2 * iRow - 1) / 2));
2918  if (dClusterDist < dClusterSize) bCheckOtherRow = kTRUE;
2919  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2920  else {
2921  // Horizontal => base = lower edge for row 0 and upper edge for row 1
2922  while (dClusterDist < dClusterSize
2923  && iRow * iNbCh / 2.0 < iMinChanInd) {
2924  iMinChanInd--;
2925  dClusterDist =
2926  TMath::Abs(poipos_local[0]
2927  - (+fChannelInfo->GetSizex()
2928  * (iMinChanInd - iChannel + (1 - 2 * iRow) / 2)));
2929  } // while right/left edge inside cluster and index not out of rpc
2930  dClusterDist = 0;
2931  while (dClusterDist < dClusterSize
2932  && iMaxChanInd < (1 + iRow) * iNbCh / 2.0 - 1) {
2933  iMaxChanInd++;
2934  dClusterDist =
2935  TMath::Abs(poipos_local[0]
2936  - (+fChannelInfo->GetSizex()
2937  * (iMaxChanInd - iChannel - (1 - 2 * iRow) / 2)));
2938  } // while left/right edge inside cluster and index not out of rpc
2939 
2940  // Test if Pad in diff row but same column as central pad has cluster overlap
2941  // if Yes => Loop from min+1 to max-1 equivalents in the opposite row
2942  dClusterDist = TMath::Abs(
2943  poipos_local[1] - (+fChannelInfo->GetSizey() * (1 - 2 * iRow) / 2));
2944  if (dClusterDist < dClusterSize) bCheckOtherRow = kTRUE;
2945  } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2946 
2947  // Loop over all channels inside the interval ]iMinChanInd;iMaxChanInd[ except central channel,
2948  // using the overlap area with cluster to get the charge and adding a Digi for all channels
2949  for (Int_t iChanInd = iMinChanInd + 1; iChanInd < iMaxChanInd;
2950  iChanInd++) {
2951  if (iChanInd == iChannel) continue;
2952 
2953  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
2954  // ... don't ask me why ... reason found in TofMC and TofGeoHandler!
2955  Int_t iCh1 = iChanInd;
2956  if (fGeoHandler->GetGeoVersion() < k14a)
2957  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
2958  CbmTofDetectorInfo xDetInfo(
2959  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
2960 
2961  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
2962 
2963  fChannelInfo = fDigiPar->GetCell(iSideChId);
2964 
2965  // Calculate the fraction of the charge in this channel
2966  Double_t dChargeSideCh =
2967  dClustCharge
2968  * ComputeClusterAreaOnChannel(
2969  iSideChId, dClusterSize, poipos_local[0], poipos_local[1]);
2970  dChargeSideCh /= dClustArea;
2971 
2972  // Fee Threshold on charge
2973  if (dChargeSideCh
2974  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iChanInd]
2975  < fDigiBdfPar->GetFeeThreshold())
2976  continue;
2977 
2978  Double_t dClustToReadout = 0.0;
2979  Double_t dPadTime = dCentralTime;
2980  // First calculate the propagation time to the center of the pad base
2981  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
2982  // Vertical => base = right/left edge
2983  dClustToReadout =
2984  TMath::Sqrt(TMath::Power(poipos_local[1], 2)
2985  + TMath::Power(poipos_local[0]
2986  - (+(1 - 2 * iRow)
2987  * fChannelInfo->GetSizex() / 2.0),
2988  2));
2989  else // Horizontal => base = bottom/upper edge
2990  dClustToReadout =
2991  TMath::Sqrt(TMath::Power(poipos_local[0], 2)
2992  + TMath::Power(poipos_local[1]
2993  - (-(1 - 2 * iRow)
2994  * fChannelInfo->GetSizey() / 2.0),
2995  2));
2996 
2997  dPadTime +=
2998  gRandom->Gaus(0.0, fdTimeResElec)
2999  + dClustToReadout / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3000 
3001  CbmTofDigi* tofDigi = new CbmTofDigi(
3002  iSM,
3003  iRpc,
3004  iChanInd,
3005  dPadTime,
3006  dChargeSideCh * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iChanInd],
3007  0,
3008  iSmType);
3009  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3010  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
3011  fStorDigi[iSmType][iSM * iNbRpc + iRpc][iChanInd].push_back(data);
3012  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iChanInd].push_back(
3013  iPntInd);
3014  LOG(debug1) << Form(
3015  "Digimatch (%d,%d,%d,%d): size %zu, MCpP %d, MCtP %d",
3016  iSmType,
3017  iSM,
3018  iRpc,
3019  iChannel,
3020  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].size(),
3021  iPntInd,
3022  iTrackID);
3023  } // for( Int_t iChanInd = iMinChanInd + 1; iChanInd < iMaxChanInd; iChanInd++ )
3024 
3025  // Loop over channels on the other row equivalent to the interval ]iMinChanInd;iMaxChanInd[
3026  if (kTRUE == bCheckOtherRow)
3027  for (Int_t iChanInd = iMinChanInd + 1 + (1 - 2 * iRow) * iNbCh / 2.0;
3028  iChanInd < iMaxChanInd + (1 - 2 * iRow) * iNbCh / 2.0;
3029  iChanInd++) {
3030  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
3031  // ... don't ask me why ... reason found in TofMC and TofGeoHandler!
3032  Int_t iCh1 = iChanInd;
3033  if (fGeoHandler->GetGeoVersion() < k14a)
3034  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
3035  CbmTofDetectorInfo xDetInfo(
3036  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
3037 
3038  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
3039 
3040  fChannelInfo = fDigiPar->GetCell(iSideChId);
3041 
3042  // Calculate the fraction of the charge in this channel
3043  Double_t dChargeSideCh =
3044  dClustCharge
3045  * ComputeClusterAreaOnChannel(
3046  iSideChId, dClusterSize, poipos_local[0], poipos_local[1]);
3047 
3048  // Fee Threshold on charge
3049  if (dChargeSideCh
3050  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iChanInd]
3051  < fDigiBdfPar->GetFeeThreshold())
3052  continue;
3053 
3054  Double_t dClustToReadout = 0.0;
3055  Double_t dPadTime = dCentralTime;
3056  // First calculate the propagation time to the center of the pad base
3057  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
3058  // Vertical => base = right/left edge
3059  dClustToReadout =
3060  TMath::Sqrt(TMath::Power(poipos_local[1], 2)
3061  + TMath::Power(poipos_local[0]
3062  - (+(1 - 2 * iRow)
3063  * fChannelInfo->GetSizex() / 2.0),
3064  2));
3065  else // Horizontal => base = bottom/upper edge
3066  dClustToReadout =
3067  TMath::Sqrt(TMath::Power(poipos_local[0], 2)
3068  + TMath::Power(poipos_local[1]
3069  - (-(1 - 2 * iRow)
3070  * fChannelInfo->GetSizey() / 2.0),
3071  2));
3072 
3073  dPadTime +=
3074  gRandom->Gaus(0.0, fdTimeResElec)
3075  + dClustToReadout / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3076 
3077  CbmTofDigi* tofDigi = new CbmTofDigi(
3078  iSM,
3079  iRpc,
3080  iChanInd,
3081  dPadTime,
3082  dChargeSideCh
3083  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iChanInd],
3084  0,
3085  iSmType);
3086  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3087  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
3088  fStorDigi[iSmType][iSM * iNbRpc + iRpc][iChanInd].push_back(data);
3089  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iChanInd].push_back(
3090  iPntInd);
3091  LOG(debug1) << Form(
3092  "Digimatch (%d,%d,%d,%d): size %zu, MCpX %d, MCtX %d",
3093  iSmType,
3094  iSM,
3095  iRpc,
3096  iChannel,
3097  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1]
3098  .size(),
3099  iPntInd,
3100  iTrackID);
3101 
3102  } // for( Int_t iChanInd = iMinChanInd + 1; iChanInd < iMaxChanInd; iChanInd++ )
3103  } // else of if( 0 == iChType)
3104 
3105  } // for (Int_t iPntInd = 0; iPntInd < nTofPoint; iPntInd++ )
3106 
3107  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
3108  // Clear the Track to channel temporary storage
3109  for (UInt_t uTrkInd = 0; uTrkInd < fvlTrckChAddr.size(); uTrkInd++) {
3110  fvlTrckChAddr[uTrkInd].clear();
3111  fvlTrckRpcAddr[uTrkInd].clear();
3112  fvlTrckRpcTime[uTrkInd].clear();
3113  }
3114  fvlTrckChAddr.clear();
3115  fvlTrckRpcAddr.clear();
3116  fvlTrckRpcTime.clear();
3117  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
3118 
3119  return kTRUE;
3120 }
3121 /************************************************************************************/
3122 // Functions for a 2D Gauss cluster approximation
3123 Bool_t CbmTofDigitize::DigitizeGaussCharge() {
3124  // Uniform distribution in ]0;x]
3125  // gRandom->Uniform(x);
3126  // Gauss distribution in of mean m, sigma s
3127  // gRandom->Gauss(m, s);
3128 
3129  CbmTofPoint* pPoint;
3130  CbmMCTrack* pMcTrk;
3131 
3132  Int_t nTofPoint = fTofPointsColl->GetEntries();
3133  Int_t nMcTracks = fMcTracksColl->GetEntries();
3134 
3135  // Prepare the temporary storing of the Track/Point/Digi info
3136  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
3137  fvlTrckChAddr.resize(nMcTracks);
3138  fvlTrckRpcAddr.resize(nMcTracks);
3139  fvlTrckRpcTime.resize(nMcTracks);
3140  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
3141  for (Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++) {
3142  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
3143  fvlTrckChAddr[iTrkInd].clear();
3144  fvlTrckRpcAddr[iTrkInd].clear();
3145  fvlTrckRpcTime[iTrkInd].clear();
3146  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
3147  } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++)
3148 
3149  if (fbMcTrkMonitor) {
3150  // Debug Info printout
3151  Int_t iNbTofTracks = 0;
3152  Int_t iNbTofTracksPrim = 0;
3153 
3154  LOG(debug1) << "CbmTofDigitize::DigitizeGaussCharge: " << nTofPoint
3155  << " points in Tof for this event with " << nMcTracks
3156  << " MC tracks ";
3157  for (Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++) {
3158  pMcTrk = (CbmMCTrack*) fMcTracksColl->At(iTrkInd);
3159  if (0 < pMcTrk->GetNPoints(ECbmModuleId::kTof)) iNbTofTracks++;
3160  if (0 < pMcTrk->GetNPoints(ECbmModuleId::kTof)
3161  && -1 == pMcTrk->GetMotherId())
3162  iNbTofTracksPrim++;
3163  } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++)
3164 
3165  //Some numbers on TOF distributions
3166  LOG(debug1) << "CbmTofDigitize::DigitizeGaussCharge: " << iNbTofTracks
3167  << " tracks in Tof ";
3168  LOG(debug1) << "CbmTofDigitize::DigitizeGaussCharge: " << iNbTofTracksPrim
3169  << " tracks in Tof from vertex";
3170  } // if( fbMcTrkMonitor )
3171 
3172  // Tof Point Info
3173  Int_t iDetId, iSmType, iSM, iRpc, iChannel, iGap;
3174  Int_t iTrackID, iChanId;
3175  TVector3 vPntPos;
3176  // Cluster Info
3177  Double_t dClusterSize;
3178  Double_t dClustCharge;
3179  // Digi
3180  // CbmTofDigi * pTofDigi; // <- Compressed digi (64 bits)
3181  // CbmTofDigi * pTofDigiExp; // <- Expanded digi
3182  // General geometry info
3183  Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
3184  Int_t iNbSm, iNbRpc, iNbCh;
3185  Int_t iChType;
3186 
3187  // Start jitter: Same contrib. for all points in same event
3188  // FIXME
3189  // Actually, event start time reconstruction should not be done in the ToF
3190  // digitizer, neither in event-based nor in time-based mode. The timestamp
3191  // of a CbmTofDigi object is not the time difference between signals in a
3192  // (hypothetical) start detector system and the MRPC wall but rather
3193  // represents the detection point in time in an MRPC measured with respect to
3194  // the start of the run.
3195  // For compatibility reasons, we continue adding a start detector jitter to
3196  // digi timestamps in event-based mode. In time-based mode, event start time
3197  // reconstruction is not considered here.
3198  Double_t dStartJitter = 0.;
3199  if (!fbTimeBasedOutput) {
3200  dStartJitter = gRandom->Gaus(0.0, fDigiBdfPar->GetStartTimeRes());
3201  }
3202 
3203  for (Int_t iPntInd = 0; iPntInd < nTofPoint; iPntInd++) {
3204  // Get a pointer to the TOF point
3205  pPoint = (CbmTofPoint*) fTofPointsColl->At(iPntInd);
3206  if (NULL == pPoint) {
3207  LOG(WARNING) << "CbmTofDigitize::DigitizeGaussCharge => Be careful: hole "
3208  "in the CbmTofPoint TClonesArray!"
3209  << endl;
3210  continue;
3211  } // if( pPoint )
3212 
3213  // Get all channel info
3214  iDetId = pPoint->GetDetectorID();
3215  iSmType = fGeoHandler->GetSMType(iDetId);
3216  iRpc = fGeoHandler->GetCounter(iDetId);
3217 
3218  iChannel = fGeoHandler->GetCell(iDetId);
3219  iGap = fGeoHandler->GetGap(iDetId);
3220  iSM = fGeoHandler->GetSModule(iDetId);
3221  iChanId = fGeoHandler->GetCellId(iDetId);
3222  iTrackID = pPoint->GetTrackID();
3223 
3224  iNbSm = fDigiBdfPar->GetNbSm(iSmType);
3225  iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
3226  iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
3227  iChType = fDigiBdfPar->GetChanType(iSmType, iRpc);
3228  Int_t iTrkId = pPoint->GetTrackID();
3229 
3230  // Catch case where the MC track was removed from the transport stack
3231  if (iTrkId < 0) {
3232  if (fbAllowPointsWithoutTrack)
3233  continue;
3234  else
3235  LOG(fatal)
3236  << "CbmTofDigitize::DigitizeDirectClusterSize => TofPoint without "
3237  "valid MC track Index, "
3238  << " track was probably cut at the transport level! "
3239  << " Pnt Idx: " << iPntInd << " Trk Idx: " << iTrkId << "\n"
3240  << "=============> To allow this kind of points and simply jump "
3241  "them, "
3242  << "call CbmTofDigitize::AllowPointsWithoutTrack() in your macro!!";
3243  } // if( iTrkId < 0 )
3244 
3245  if (fbMcTrkMonitor) {
3246  // Get pointer to the MC-Track info
3247  pMcTrk = (CbmMCTrack*) fMcTracksColl->At(iTrkId);
3248 
3249  // Keep only primaries
3250  if (kTRUE == fDigiBdfPar->UseOnlyPrimaries())
3251  if (-1 != pMcTrk->GetMotherId()) continue;
3252  } // if( fbMcTrkMonitor )
3253 
3254  if (iSmType < 0. || iNbSmTypes < iSmType || iSM < 0. || iNbSm < iSM
3255  || iRpc < 0. || iNbRpc < iRpc || iChannel < 0. || iNbCh < iChannel) {
3256  LOG(error) << Form("CbmTofDigitize => det ID 0x%08x", iDetId)
3257  << " SMType: " << iSmType;
3258  LOG(error) << " SModule: " << iSM << " of " << iNbSm + 1;
3259  LOG(error) << " Module: " << iRpc << " of " << iNbRpc + 1;
3260  LOG(error) << " Gap: " << iGap;
3261  LOG(error) << " Cell: " << iChannel << " of " << iNbCh + 1;
3262  continue;
3263  } // if error on channel data
3264 
3265  // Check if there was already a Digi from the same track created
3266  // with this channel as main channel
3267  ULong64_t uAddr =
3268  CbmTofAddress::GetUniqueAddress(iSM, iRpc, iChannel, 0, iSmType);
3269  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
3270  Bool_t bFoundIt = kFALSE;
3271  // Check if this track ID is bigger than size of McTracks TClonesArray (maybe happen with PSD cut!!)
3272  // If yes, resize the array to reach the appropriate number of fields
3273  // Cast of TrackId is safe as we already check for "< 0" case
3274  // TODO: Is this check still required when Id < 0 are jumped?
3275  if (fvlTrckChAddr.size() <= static_cast<UInt_t>(iTrkId))
3276  fvlTrckChAddr.resize(iTrkId + 1);
3277 
3278  for (UInt_t uTrkMainCh = 0; uTrkMainCh < fvlTrckChAddr[iTrkId].size();
3279  uTrkMainCh++)
3280  if (uAddr == fvlTrckChAddr[iTrkId][uTrkMainCh]) {
3281  bFoundIt = kTRUE;
3282  break;
3283  }
3284  // If it is the case, no need to redo the full stuff for this gap
3285  if (kTRUE == bFoundIt) continue;
3286  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
3287 
3288  // Probability that the RPC detect the hit
3289  // if( fDigiBdfPar->GetEfficiency(iSmType) < gRandom->Uniform(1) )
3290  // continue;
3291 
3292  // Probability that the gap detect the hit
3293  // For now use a simple gap treatment (cf CbmTofDigiBdfPar):
3294  // - a single fired gap fires the channel
3295  // - all gaps have exactly equivalent efficiency/firing probability
3296  // - independent gap firing (no charge accumulation or such)
3297  // The probability for a hit not to fire at all is then
3298  // (1-p)^NGaps with p = gap efficiency and Ngaps the number of gaps in this RPC
3299  // This results in the relation: p = 1 - (1 - P)^(1/Ngaps)
3300  // with P = RPC efficiency from beamtime
3301  if (fDigiBdfPar->GetGapEfficiency(iSmType, iRpc) < gRandom->Uniform(1))
3302  continue;
3303 
3304  // Save the address in vector so that this channel is not redone later for the
3305  // eventual other gaps TofPoint
3306  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk())
3307  fvlTrckChAddr[iTrkId].push_back(uAddr);
3308 
3309  // Get TofPoint Position
3310  pPoint->Position(vPntPos);
3311  fChannelInfo = fDigiPar->GetCell(iChanId);
3312  TGeoNode* fNode = // prepare global->local trafo
3313  gGeoManager->FindNode(
3314  fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
3315  LOG(debug1) << Form(" TofDigitizerBDF:: (%3d,%3d,%3d,%3d) - node at "
3316  "(%6.1f,%6.1f,%6.1f) : 0x%p Pos(%6.1f,%6.1f,%6.1f)",
3317  iSmType,
3318  iSM,
3319  iRpc,
3320  iChannel,
3321  fChannelInfo->GetX(),
3322  fChannelInfo->GetY(),
3323  fChannelInfo->GetZ(),
3324  fNode,
3325  vPntPos.X(),
3326  vPntPos.Y(),
3327  vPntPos.Z());
3328  gGeoManager->GetCurrentNode();
3329  gGeoManager->GetCurrentMatrix();
3330  Double_t poipos[3] = {vPntPos.X(), vPntPos.Y(), vPntPos.Z()};
3331  Double_t poipos_local[3];
3332  gGeoManager->MasterToLocal(poipos, poipos_local);
3333 
3334  // Generate a Cluster radius in [cm]
3335  dClusterSize = GenerateClusterRadius(iSmType, iRpc);
3336  while (dClusterSize < 0.0001)
3337  // Should not happen without error message
3338  // But Landau can give really small values
3339  dClusterSize = GenerateClusterRadius(iSmType, iRpc);
3340 
3341  // Generate Cluster charge as ToT[ps]
3342  dClustCharge = fh1ClusterTotProb[iSmType]->GetBinContent(
3343  fh1ClusterTotProb[iSmType]->FindBin(gRandom->Uniform(1)));
3344 
3345  // Assume the width (sigma) of the gaussian in both direction is dClusterSize/3
3346  // => 99% of the charge is in "dClusterSize"
3347  TF2* f2ChargeDist =
3348  new TF2("ChargeDist",
3349  "[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])",
3350  -5.0 * dClusterSize,
3351  5.0 * dClusterSize,
3352  -5.0 * dClusterSize,
3353  5.0 * dClusterSize);
3354  f2ChargeDist->SetParameters(
3355  dClustCharge / (TMath::Sqrt(2.0 * TMath::Pi()) * dClusterSize / 3.0),
3356  poipos_local[0],
3357  dClusterSize / 3.0,
3358  poipos_local[1],
3359  dClusterSize / 3.0);
3360 
3361  // Calculate the fraction of the charge in central channel
3362  Double_t dChargeCentral = f2ChargeDist->Integral(
3363  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
3364  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
3365  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
3366  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
3367 
3368  // Calculate the time for the central channel
3369  Double_t dCentralTime;
3370  // Check if there was already a Digi from the same track created in this RPC
3371  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
3372  ULong64_t uRpcAddr =
3373  CbmTofAddress::GetUniqueAddress(iSM, iRpc, 0, 0, iSmType);
3374  Bool_t bFoundIt = kFALSE;
3375  UInt_t uTrkRpcPair = 0;
3376  for (uTrkRpcPair = 0; uTrkRpcPair < fvlTrckRpcAddr[iTrkId].size();
3377  uTrkRpcPair++)
3378  if (uRpcAddr == fvlTrckRpcAddr[iTrkId][uTrkRpcPair]) {
3379  bFoundIt = kTRUE;
3380  break;
3381  }
3382  // If it is the case, we should reuse the timing already assigned to this track
3383  if (kTRUE == bFoundIt)
3384  dCentralTime = fvlTrckRpcTime[iTrkId][uTrkRpcPair];
3385  else {
3386  dCentralTime =
3387  pPoint->GetTime()
3388  + gRandom->Gaus(0.0, fDigiBdfPar->GetResolution(iSmType))
3389  + dStartJitter; // Same contrib. for all points in same event
3390  fvlTrckRpcAddr[iTrkId].push_back(uRpcAddr);
3391  fvlTrckRpcTime[iTrkId].push_back(dCentralTime);
3392  } // No Digi yet in this RPC for this Track
3393  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
3394  else
3395  dCentralTime =
3396  pPoint->GetTime()
3397  + gRandom->Gaus(0.0, fDigiBdfPar->GetResolution(iSmType))
3398  + dStartJitter; // Same contrib. for all points in same event
3399 
3400  // Calculate propagation time(s) to the readout point(s)
3401  if (0 == iChType) {
3402  // Strips (readout on both side)
3403  // Assume that the bottom/left strip have lower channel index
3404  // E.g: ... | i | i+1 | ...
3405  // or ... y
3406  // ----- ^
3407  // i+1 |
3408  // ----- --> x
3409  // i
3410  // -----
3411  // ...
3412 
3413  Double_t dTimeA = dCentralTime;
3414  Double_t dTimeB = dCentralTime;
3415  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
3416  // Horizontal channel: A = Right and B = Left of the channel
3417  dTimeA +=
3418  gRandom->Gaus(0.0, fdTimeResElec)
3419 #ifdef FULL_PROPAGATION_TIME
3420  + TMath::Abs(poipos_local[0] - (+fChannelInfo->GetSizex() / 2.0))
3421 #else
3422  + (poipos_local[0])
3423 #endif
3424  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3425  dTimeB +=
3426  gRandom->Gaus(0.0, fdTimeResElec)
3427 #ifdef FULL_PROPAGATION_TIME
3428  + TMath::Abs(poipos_local[0] - (-fChannelInfo->GetSizex() / 2.0))
3429 #else
3430  - (poipos_local[0])
3431 #endif
3432  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3433  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
3434  else {
3435  // Vertical channel: A = Top and B = Bottom of the channel
3436  dTimeA +=
3437  gRandom->Gaus(0.0, fdTimeResElec)
3438 #ifdef FULL_PROPAGATION_TIME
3439  + TMath::Abs(poipos_local[1] - (+fChannelInfo->GetSizey() / 2.0))
3440 #else
3441  + (poipos_local[1])
3442 #endif
3443  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3444  dTimeB +=
3445  gRandom->Gaus(0.0, fdTimeResElec)
3446 #ifdef FULL_PROPAGATION_TIME
3447  + TMath::Abs(poipos_local[1] - (-fChannelInfo->GetSizey() / 2.0))
3448 #else
3449  - (poipos_local[1])
3450 #endif
3451  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3452  } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
3453 
3454  CbmTofDigi* tofDigi = new CbmTofDigi(
3455  iSM,
3456  iRpc,
3457  iChannel,
3458  dTimeA,
3459  dChargeCentral
3460  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1] / 2.0,
3461  1,
3462  iSmType);
3463  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3464  std::pair<CbmTofDigi*, CbmMatch*> data1(tofDigi, nullptr);
3465  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].push_back(
3466  data1);
3467  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].push_back(
3468  iPntInd);
3469  LOG(debug1) << Form(
3470  "Digimatch (%d,%d,%d,%d): size %zu, val1 %d, MCt %d",
3471  iSmType,
3472  iSM,
3473  iRpc,
3474  iChannel,
3475  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel + 1].size(),
3476  iPntInd,
3477  iTrackID);
3478 
3479  tofDigi = new CbmTofDigi(
3480  iSM,
3481  iRpc,
3482  iChannel,
3483  dTimeB,
3484  dChargeCentral
3485  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iChannel] / 2.0,
3486  0,
3487  iSmType);
3488  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3489  std::pair<CbmTofDigi*, CbmMatch*> data2(tofDigi, nullptr);
3490  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iChannel].push_back(data2);
3491  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel].push_back(
3492  iPntInd);
3493  LOG(debug1) << Form(
3494  "Digimatch (%d,%d,%d,%d): size %zu, val2 %d, MCt %d",
3495  iSmType,
3496  iSM,
3497  iRpc,
3498  iChannel,
3499  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iChannel].size(),
3500  iPntInd,
3501  iTrackID);
3502  } // if( 0 == iChType)
3503  else {
3504  // Pad (Single side readout)
3505  // Assume that the bottom/left pads have lower channel index
3506  // E.g: for a 2*4 pads counter, pads 0-3 are the bottom/left ones and 4-7 the
3507  // top/right one with reversed numbering:
3508  // -----------------
3509  // row 1 | 7 | 6 | 5 | 4 | y
3510  // ----------------- ^
3511  // row 0 | 0 | 1 | 2 | 3 | |
3512  // ----------------- --> x
3513  // or ---------
3514  // | 4 | 3 |
3515  // ---------
3516  // | 5 | 2 |
3517  // ---------
3518  // | 6 | 1 |
3519  // ---------
3520  // | 7 | 0 |
3521  // ---------
3522  // row 1 0
3523  // Also assume that the readout happens at the middle of the readout edge
3524 
3525  Double_t dClustToReadout = 0.0;
3526  Double_t dPadTime = dCentralTime;
3527  // First calculate the propagation time to the center of the pad base
3528  if (iChannel < iNbCh / 2.0) {
3529  // First row
3530  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
3531  // Vertical => base = right edge
3532  dClustToReadout = TMath::Sqrt(
3533  TMath::Power(poipos_local[1], 2)
3534  + TMath::Power(poipos_local[0] - (fChannelInfo->GetSizex() / 2.0),
3535  2));
3536  else // Horizontal => base = bottom edge
3537  dClustToReadout = TMath::Sqrt(
3538  TMath::Power(poipos_local[0], 2)
3539  + TMath::Power(poipos_local[1] - (-fChannelInfo->GetSizey() / 2.0),
3540  2));
3541  } // if( iChannel < iNbCh/2.0 )
3542  else {
3543  // Second row
3544  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
3545  // Vertical => base = left edge
3546  dClustToReadout = TMath::Sqrt(
3547  TMath::Power(poipos_local[1], 2)
3548  + TMath::Power(poipos_local[0] - (-fChannelInfo->GetSizex() / 2.0),
3549  2));
3550  else // Horizontal => base = upper edge
3551  dClustToReadout = TMath::Sqrt(
3552  TMath::Power(poipos_local[0], 2)
3553  + TMath::Power(poipos_local[1] - (+fChannelInfo->GetSizey() / 2.0),
3554  2));
3555  } // else of if( iChannel < iNbCh/2.0 )
3556 
3557  dPadTime += gRandom->Gaus(0.0, fdTimeResElec)
3558  + dClustToReadout / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3559 
3560  // TODO: Check on fee threshold ?
3561  CbmTofDigi* tofDigi = new CbmTofDigi(
3562  iSM,
3563  iRpc,
3564  iChannel,
3565  dPadTime,
3566  dChargeCentral * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iChannel],
3567  0,
3568  iSmType);
3569  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3570  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
3571  fStorDigi[iSmType][iSM * iNbRpc + iRpc][iChannel].push_back(data);
3572  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iChannel].push_back(iPntInd);
3573  LOG(debug1) << Form(
3574  "Digimatch (%d,%d,%d,%d): size %zu, val3 %d, MCt %d",
3575  iSmType,
3576  iSM,
3577  iRpc,
3578  iChannel,
3579  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iChannel].size(),
3580  iPntInd,
3581  iTrackID);
3582  } // else of if( 0 == iChType)
3583 
3584  // Loop over neighboring channel and check if they are above threshold
3585  if (0 == iChType) {
3586  // Strips (readout on both side)
3587  // Loop in decreasing channel index until a channel with charge under threshold is found
3588  if (0 < iChannel) {
3589  Int_t iSideChInd = iChannel - 1;
3590 
3591  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
3592  // ... don't ask me why ... reason found in TofMC and TofGeoHandler!
3593  Int_t iCh1 = iSideChInd;
3594  if (fGeoHandler->GetGeoVersion() < k14a)
3595  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
3596  CbmTofDetectorInfo xDetInfo(
3597  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
3598 
3599  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
3600 
3601  fChannelInfo = fDigiPar->GetCell(iSideChId);
3602 
3603  // Calculate the fraction of the charge in this channel
3604  Double_t dChargeSideCh = f2ChargeDist->Integral(
3605  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
3606  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
3607  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
3608  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
3609 
3610  while (fDigiBdfPar->GetFeeThreshold() <= dChargeSideCh
3611  && 0 <= iSideChInd) {
3612  // Strips = readout on both side
3613  Double_t dTimeA = dCentralTime;
3614  Double_t dTimeB = dCentralTime;
3615  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
3616  // Horizontal channel: A = Right and B = Left of the channel
3617  dTimeA +=
3618  gRandom->Gaus(0.0, fdTimeResElec)
3619 #ifdef FULL_PROPAGATION_TIME
3620  + TMath::Abs(poipos_local[0] - (+fChannelInfo->GetSizex() / 2.0))
3621 #else
3622  + (poipos_local[0])
3623 #endif
3624  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3625  dTimeB +=
3626  gRandom->Gaus(0.0, fdTimeResElec)
3627 #ifdef FULL_PROPAGATION_TIME
3628  + TMath::Abs(poipos_local[0] - (-fChannelInfo->GetSizex() / 2.0))
3629 #else
3630  - (poipos_local[0])
3631 #endif
3632  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3633  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
3634  else {
3635  // Vertical channel: A = Top and B = Bottom of the channel
3636  dTimeA +=
3637  gRandom->Gaus(0.0, fdTimeResElec)
3638 #ifdef FULL_PROPAGATION_TIME
3639  + TMath::Abs(poipos_local[1] - (+fChannelInfo->GetSizey() / 2.0))
3640 #else
3641  + (poipos_local[1])
3642 #endif
3643  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3644  dTimeB +=
3645  gRandom->Gaus(0.0, fdTimeResElec)
3646 #ifdef FULL_PROPAGATION_TIME
3647  + TMath::Abs(poipos_local[1] - (-fChannelInfo->GetSizey() / 2.0))
3648 #else
3649  - (poipos_local[1])
3650 #endif
3651  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3652  } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
3653 
3654  CbmTofDigi* tofDigi = new CbmTofDigi(
3655  iSM,
3656  iRpc,
3657  iSideChInd,
3658  dTimeA,
3659  dChargeSideCh
3660  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd + 1]
3661  / 2.0,
3662  1,
3663  iSmType);
3664  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3665  std::pair<CbmTofDigi*, CbmMatch*> data1(tofDigi, nullptr);
3666  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd + 1].push_back(
3667  data1);
3668  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd + 1]
3669  .push_back(iPntInd);
3670  tofDigi = new CbmTofDigi(
3671  iSM,
3672  iRpc,
3673  iSideChInd,
3674  dTimeB,
3675  dChargeSideCh
3676  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd]
3677  / 2.0,
3678  0,
3679  iSmType);
3680  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3681  std::pair<CbmTofDigi*, CbmMatch*> data2(tofDigi, nullptr);
3682  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd].push_back(
3683  data2);
3684  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd]
3685  .push_back(iPntInd);
3686  LOG(debug1) << Form(
3687  "Digimatch (%d,%d,%d,%d): size %zu, val4 %d, MCt %d",
3688  iSmType,
3689  iSM,
3690  iRpc,
3691  iChannel,
3692  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd].size(),
3693  iPntInd,
3694  iTrackID);
3695 
3696 
3697  // Check next
3698  iSideChInd--;
3699 
3700  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
3701  // ... don't ask me why ...
3702  // FIXME: problem around here => reason found in TofMC and TofGeoHandler!
3703  Int_t iCh2 = iSideChInd;
3704  if (fGeoHandler->GetGeoVersion() < k14a)
3705  iCh2 = iCh2 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
3706  xDetInfo.fCell = iCh2;
3707 
3708  iSideChId = fTofId->SetDetectorInfo(xDetInfo);
3709 
3710  fChannelInfo = fDigiPar->GetCell(iSideChId);
3711 
3712  // Calculate the fraction of the charge in this channel
3713  dChargeSideCh = f2ChargeDist->Integral(
3714  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
3715  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
3716  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
3717  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
3718  } // while signal and channel index ok
3719  } // if( 0 < iChannel)
3720  // Loop in increasing channel index until a channel with charge under threshold is found
3721  if (iChannel < iNbCh - 1) {
3722  Int_t iSideChInd = iChannel + 1;
3723 
3724  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
3725  // ... don't ask me why ... reason found in TofMC and TofGeoHandler!
3726  Int_t iCh1 = iSideChInd;
3727  if (fGeoHandler->GetGeoVersion() < k14a)
3728  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
3729  CbmTofDetectorInfo xDetInfo(
3730  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
3731 
3732  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
3733 
3734  fChannelInfo = fDigiPar->GetCell(iSideChId);
3735 
3736  // Calculate the fraction of the charge in this channel
3737  Double_t dChargeSideCh = f2ChargeDist->Integral(
3738  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
3739  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
3740  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
3741  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
3742 
3743  while (fDigiBdfPar->GetFeeThreshold() < dChargeSideCh
3744  && iSideChInd < iNbCh) {
3745  // Strips = readout on both side
3746  Double_t dTimeA = dCentralTime;
3747  Double_t dTimeB = dCentralTime;
3748  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
3749  // Horizontal channel: A = Right and B = Left of the channel
3750  dTimeA +=
3751  gRandom->Gaus(0.0, fdTimeResElec)
3752 #ifdef FULL_PROPAGATION_TIME
3753  + TMath::Abs(poipos_local[0] - (+fChannelInfo->GetSizex() / 2.0))
3754 #else
3755  + (poipos_local[0])
3756 #endif
3757  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3758  dTimeB +=
3759  gRandom->Gaus(0.0, fdTimeResElec)
3760 #ifdef FULL_PROPAGATION_TIME
3761  + TMath::Abs(poipos_local[0] - (-fChannelInfo->GetSizex() / 2.0))
3762 #else
3763  - (poipos_local[0])
3764 #endif
3765  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3766  } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
3767  else {
3768  // Vertical channel: A = Top and B = Bottom of the channel
3769  dTimeA +=
3770  gRandom->Gaus(0.0, fdTimeResElec)
3771 #ifdef FULL_PROPAGATION_TIME
3772  + TMath::Abs(poipos_local[1] - (+fChannelInfo->GetSizey() / 2.0))
3773 #else
3774  + (poipos_local[1])
3775 #endif
3776  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3777  dTimeB +=
3778  gRandom->Gaus(0.0, fdTimeResElec)
3779 #ifdef FULL_PROPAGATION_TIME
3780  + TMath::Abs(poipos_local[1] - (-fChannelInfo->GetSizey() / 2.0))
3781 #else
3782  - (poipos_local[1])
3783 #endif
3784  / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3785  } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
3786 
3787  CbmTofDigi* tofDigi = new CbmTofDigi(
3788  iSM,
3789  iRpc,
3790  iSideChInd,
3791  dTimeA,
3792  dChargeSideCh
3793  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd + 1]
3794  / 2.0,
3795  1,
3796  iSmType);
3797  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3798  std::pair<CbmTofDigi*, CbmMatch*> data1(tofDigi, nullptr);
3799  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd + 1].push_back(
3800  data1);
3801  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd + 1]
3802  .push_back(iPntInd);
3803  LOG(debug1) << Form(
3804  "Digimatch (%d,%d,%d,%d): size %zu, val5 %d, MCt %d",
3805  iSmType,
3806  iSM,
3807  iRpc,
3808  iChannel,
3809  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd + 1]
3810  .size(),
3811  iPntInd,
3812  iTrackID);
3813 
3814  tofDigi = new CbmTofDigi(
3815  iSM,
3816  iRpc,
3817  iSideChInd,
3818  dTimeB,
3819  dChargeSideCh
3820  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd]
3821  / 2.0,
3822  0,
3823  iSmType);
3824  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3825  std::pair<CbmTofDigi*, CbmMatch*> data2(tofDigi, nullptr);
3826  fStorDigi[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd].push_back(
3827  data2);
3828  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd]
3829  .push_back(iPntInd);
3830  LOG(debug1) << Form(
3831  "Digimatch (%d,%d,%d,%d): size %zu, val6 %d, MCt %d",
3832  iSmType,
3833  iSM,
3834  iRpc,
3835  iChannel,
3836  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd].size(),
3837  iPntInd,
3838  iTrackID);
3839 
3840  // Check next
3841  iSideChInd++;
3842 
3843  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
3844  // ... don't ask me why ...
3845  xDetInfo.fCell = iSideChInd + 1;
3846 
3847  iSideChId = fTofId->SetDetectorInfo(xDetInfo);
3848 
3849  fChannelInfo = fDigiPar->GetCell(iSideChId);
3850 
3851  // Calculate the fraction of the charge in this channel
3852  dChargeSideCh = f2ChargeDist->Integral(
3853  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
3854  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
3855  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
3856  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
3857  } // while signal and channel index ok
3858  } // if( 0 < iChannel)
3859  } // if( 0 == iChType)
3860  else {
3861  // Pad (Single side readout)
3862  // Assume that the bottom/left pads have lower channel index
3863  // E.g: for a 2*4 pads counter, pads 0-3 are the bottom/left ones and 4-7 the
3864  // top/right one with reversed numbering
3865  // -----------------
3866  // row 1 | 7 | 6 | 5 | 4 | y
3867  // ----------------- ^
3868  // row 0 | 0 | 1 | 2 | 3 | |
3869  // ----------------- --> x
3870  // or ---------
3871  // | 4 | 3 |
3872  // ---------
3873  // | 5 | 2 |
3874  // ---------
3875  // | 6 | 1 |
3876  // ---------
3877  // | 7 | 0 |
3878  // ---------
3879  // row 1 0
3880  // Loop in decreasing channel index in same row, find min = 1st 0 channel
3881  // Loop in increasing channel index in same row, find max = 1st 0 channel
3882  Int_t iMinChanInd = iChannel;
3883  Int_t iMaxChanInd = iChannel;
3884 
3885  // Double_t dClusterDist = 0;// -> Comment to remove warning because set but never used
3886  Int_t iRow;
3887  // Bool_t bCheckOtherRow = kFALSE; // -> Comment to remove warning because set but never used
3888  if (iChannel < iNbCh / 2.0)
3889  iRow = 0;
3890  else
3891  iRow = 1;
3892 
3893  // Loop in decreasing channel index until a channel with charge under threshold is found
3894  if (iRow * iNbCh / 2.0 < iChannel) {
3895  Int_t iSideChInd = iChannel - 1;
3896 
3897  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
3898  // ... don't ask me why ... reason found in TofMC and TofGeoHandler!
3899  Int_t iCh1 = iSideChInd;
3900  if (fGeoHandler->GetGeoVersion() < k14a)
3901  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
3902  CbmTofDetectorInfo xDetInfo(
3903  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
3904 
3905  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
3906 
3907  fChannelInfo = fDigiPar->GetCell(iSideChId);
3908 
3909  // Calculate the fraction of the charge in this channel
3910  Double_t dChargeSideCh = f2ChargeDist->Integral(
3911  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
3912  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
3913  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
3914  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
3915 
3916  while (fDigiBdfPar->GetFeeThreshold() < dChargeSideCh
3917  && iRow * iNbCh / 2.0 <= iSideChInd) {
3918  iMinChanInd = iSideChInd;
3919 
3920  Double_t dClustToReadout = 0.0;
3921  Double_t dPadTime = dCentralTime;
3922  // First calculate the propagation time to the center of the pad base
3923  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
3924  // Vertical => base = right/left edge
3925  dClustToReadout =
3926  TMath::Sqrt(TMath::Power(poipos_local[1], 2)
3927  + TMath::Power(poipos_local[0]
3928  - ((1 - 2 * iRow)
3929  * fChannelInfo->GetSizex() / 2.0),
3930  2));
3931  else // Horizontal => base = bottom/upper edge
3932  dClustToReadout =
3933  TMath::Sqrt(TMath::Power(poipos_local[0], 2)
3934  + TMath::Power(poipos_local[1]
3935  - (-(1 - 2 * iRow)
3936  * fChannelInfo->GetSizey() / 2.0),
3937  2));
3938 
3939  dPadTime +=
3940  gRandom->Gaus(0.0, fdTimeResElec)
3941  + dClustToReadout / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
3942 
3943  CbmTofDigi* tofDigi = new CbmTofDigi(
3944  iSM,
3945  iRpc,
3946  iSideChInd,
3947  dPadTime,
3948  dChargeSideCh
3949  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iSideChInd],
3950  0,
3951  iSmType);
3952  // tofDigi->GetMatch()->AddLink(1., iPntInd);
3953  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
3954  fStorDigi[iSmType][iSM * iNbRpc + iRpc][iSideChInd].push_back(data);
3955  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iSideChInd].push_back(
3956  iPntInd);
3957  LOG(debug1) << Form(
3958  "Digimatch (%d,%d,%d,%d): size %zu, val7 %d, MCt %d",
3959  iSmType,
3960  iSM,
3961  iRpc,
3962  iChannel,
3963  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][2 * iSideChInd].size(),
3964  iPntInd,
3965  iTrackID);
3966 
3967 
3968  // Check next
3969  iSideChInd--;
3970 
3971  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
3972  // ... don't ask me why ...
3973  xDetInfo.fCell = iSideChInd + 1;
3974 
3975  iSideChId = fTofId->SetDetectorInfo(xDetInfo);
3976 
3977  fChannelInfo = fDigiPar->GetCell(iSideChId);
3978 
3979  // Calculate the fraction of the charge in this channel
3980  dChargeSideCh = f2ChargeDist->Integral(
3981  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
3982  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
3983  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
3984  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
3985  } // while channel has good charge and is not out
3986  } // if( iRow*iNbCh/2.0 < iChannel )
3987  // Loop in increasing channel index until a channel with charge under threshold is found
3988  if (iChannel < (1 + iRow) * iNbCh / 2.0 - 1) {
3989  Int_t iSideChInd = iChannel + 1;
3990 
3991  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
3992  // ... don't ask me why ... reason found in TofMC and TofGeoHandler!
3993  Int_t iCh1 = iSideChInd;
3994  if (fGeoHandler->GetGeoVersion() < k14a)
3995  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
3996  CbmTofDetectorInfo xDetInfo(
3997  ECbmModuleId::kTof, iSmType, iSM, iRpc, 0, iCh1);
3998 
3999  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
4000 
4001  fChannelInfo = fDigiPar->GetCell(iSideChId);
4002 
4003  // Calculate the fraction of the charge in this channel
4004  Double_t dChargeSideCh = f2ChargeDist->Integral(
4005  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
4006  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
4007  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
4008  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
4009 
4010  while (fDigiBdfPar->GetFeeThreshold() < dChargeSideCh
4011  && iSideChInd < (1 + iRow) * iNbCh / 2.0) {
4012  iMaxChanInd = iSideChInd;
4013 
4014  Double_t dClustToReadout = 0.0;
4015  Double_t dPadTime = dCentralTime;
4016  // First calculate the propagation time to the center of the pad base
4017  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
4018  // Vertical => base = right/left edge
4019  dClustToReadout =
4020  TMath::Sqrt(TMath::Power(poipos_local[1], 2)
4021  + TMath::Power(poipos_local[0]
4022  - (+(1 - 2 * iRow)
4023  * fChannelInfo->GetSizex() / 2.0),
4024  2));
4025  else // Horizontal => base = bottom/upper edge
4026  dClustToReadout =
4027  TMath::Sqrt(TMath::Power(poipos_local[0], 2)
4028  + TMath::Power(poipos_local[1]
4029  - (-(1 - 2 * iRow)
4030  * fChannelInfo->GetSizey() / 2.0),
4031  2));
4032 
4033  dPadTime +=
4034  gRandom->Gaus(0.0, fdTimeResElec)
4035  + dClustToReadout / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
4036 
4037  CbmTofDigi* tofDigi = new CbmTofDigi(
4038  iSM,
4039  iRpc,
4040  iSideChInd,
4041  dPadTime,
4042  dChargeSideCh
4043  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iSideChInd],
4044  0,
4045  iSmType);
4046  // tofDigi->GetMatch()->AddLink(1., iPntInd);
4047  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
4048  fStorDigi[iSmType][iSM * iNbRpc + iRpc][iSideChInd].push_back(data);
4049  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iSideChInd].push_back(
4050  iPntInd);
4051  LOG(debug1) << Form(
4052  "Digimatch (%d,%d,%d,%d): size %zu, val8 %d, MCt %d",
4053  iSmType,
4054  iSM,
4055  iRpc,
4056  iChannel,
4057  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iSideChInd].size(),
4058  iPntInd,
4059  iTrackID);
4060 
4061  // Check next
4062  iSideChInd++;
4063 
4064  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
4065  // ... don't ask me why ...
4066  xDetInfo.fCell = iSideChInd + 1;
4067 
4068  iSideChId = fTofId->SetDetectorInfo(xDetInfo);
4069 
4070  fChannelInfo = fDigiPar->GetCell(iSideChId);
4071 
4072  // Calculate the fraction of the charge in this channel
4073  dChargeSideCh = f2ChargeDist->Integral(
4074  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
4075  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
4076  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
4077  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
4078  } // while channel has good charge and is not out
4079  } // if( iChannel < (1+iRow)*iNbCh/2.0 - 1)
4080 
4081  // Test if Pad in diff row but same column as central pad has enough charge
4082  // if Yes => Loop from min to max equivalents in the opposite row
4083 
4084  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
4085  // ... don't ask me why ... reason found in TofMC and TofGeoHandler!
4086  Int_t iCh1 = iChannel;
4087  if (fGeoHandler->GetGeoVersion() < k14a)
4088  iCh1 = iCh1 + 1; //FIXME: Reason found in TofMC and TofGeoHandler
4089  CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof,
4090  iSmType,
4091  iSM,
4092  iRpc,
4093  0,
4094  iCh1 + (1 - 2 * iRow) * iNbCh / 2.0);
4095  Int_t iSideChId = fTofId->SetDetectorInfo(xDetInfo);
4096  fChannelInfo = fDigiPar->GetCell(iSideChId);
4097 
4098  // Calculate the fraction of the charge in this channel
4099  Double_t dChargeSideCh = f2ChargeDist->Integral(
4100  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
4101  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
4102  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
4103  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
4104 
4105  if (fDigiBdfPar->GetFeeThreshold() < dChargeSideCh) {
4106  // bCheckOtherRow = kTRUE; // -> Comment to remove warning because set but never used
4107 
4108  for (Int_t iChanInd = iMinChanInd + (1 - 2 * iRow) * iNbCh / 2.0;
4109  iChanInd <= iMaxChanInd + (1 - 2 * iRow) * iNbCh / 2.0;
4110  iChanInd++) {
4111  // Channel index in this UID is in [1,nbCh] instead of [0, nbCh[
4112  // ... don't ask me why ...
4113  xDetInfo.fCell = iChanInd + 1;
4114  iSideChId = fTofId->SetDetectorInfo(xDetInfo);
4115  fChannelInfo = fDigiPar->GetCell(iSideChId);
4116 
4117  // Calculate the fraction of the charge in this channel
4118  dChargeSideCh = f2ChargeDist->Integral(
4119  fChannelInfo->GetX() - fChannelInfo->GetSizex() / 2.0,
4120  fChannelInfo->GetX() + fChannelInfo->GetSizex() / 2.0,
4121  fChannelInfo->GetY() - fChannelInfo->GetSizey() / 2.0,
4122  fChannelInfo->GetY() + fChannelInfo->GetSizey() / 2.0);
4123 
4124  if (fDigiBdfPar->GetFeeThreshold() < dChargeSideCh) {
4125  Double_t dClustToReadout = 0.0;
4126  Double_t dPadTime = dCentralTime;
4127  // First calculate the propagation time to the center of the pad base
4128  if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc))
4129  // Vertical => base = right/left edge
4130  dClustToReadout = TMath::Sqrt(
4131  TMath::Power(poipos_local[1], 2)
4132  + TMath::Power(
4133  poipos_local[0]
4134  - (+(1 - 2 * iRow) * fChannelInfo->GetSizex() / 2.0),
4135  2));
4136  else // Horizontal => base = bottom/upper edge
4137  dClustToReadout = TMath::Sqrt(
4138  TMath::Power(poipos_local[0], 2)
4139  + TMath::Power(
4140  poipos_local[1]
4141  - (-(1 - 2 * iRow) * fChannelInfo->GetSizey() / 2.0),
4142  2));
4143 
4144  dPadTime +=
4145  gRandom->Gaus(0.0, fdTimeResElec)
4146  + dClustToReadout / fvdSignalVelocityRpc[iSmType][iSM][iRpc];
4147 
4148  CbmTofDigi* tofDigi = new CbmTofDigi(
4149  iSM,
4150  iRpc,
4151  iChanInd,
4152  dPadTime,
4153  dChargeSideCh
4154  * fdChannelGain[iSmType][iSM * iNbRpc + iRpc][iChanInd],
4155  0,
4156  iSmType);
4157  // tofDigi->GetMatch()->AddLink(1., iPntInd);
4158  std::pair<CbmTofDigi*, CbmMatch*> data(tofDigi, nullptr);
4159  fStorDigi[iSmType][iSM * iNbRpc + iRpc][iChanInd].push_back(data);
4160  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iChanInd].push_back(
4161  iPntInd);
4162  LOG(debug1) << Form(
4163  "Digimatch (%d,%d,%d,%d): size %zu, val9 %d, MCt %d",
4164  iSmType,
4165  iSM,
4166  iRpc,
4167  iChannel,
4168  fStorDigiMatch[iSmType][iSM * iNbRpc + iRpc][iChanInd].size(),
4169  iPntInd,
4170  iTrackID);
4171 
4172  } // if( fDigiBdfPar->GetFeeThreshold() < dChargeSideCh )
4173  } // for channels on other row where same row had signal
4174  } // if( fDigiBdfPar->GetFeeThreshold() < dChargeSideCh )
4175  } // else of if( 0 == iChType)
4176 
4177  delete f2ChargeDist;
4178  } // for (Int_t iPntInd = 0; iPntInd < nTofPoint; iPntInd++ )
4179 
4180  if (kTRUE == fDigiBdfPar->UseOneGapPerTrk()) {
4181  // Clear the Track to channel temporary storage
4182  for (UInt_t uTrkInd = 0; uTrkInd < fvlTrckChAddr.size(); uTrkInd++) {
4183  fvlTrckChAddr[uTrkInd].clear();
4184  fvlTrckRpcAddr[uTrkInd].clear();
4185  fvlTrckRpcTime[uTrkInd].clear();
4186  }
4187  fvlTrckChAddr.clear();
4188  fvlTrckRpcAddr.clear();
4189  fvlTrckRpcTime.clear();
4190  } // if( kTRUE == fDigiBdfPar->UseOneGapPerTrk() )
4191 
4192  return kTRUE;
4193 }
4194 /************************************************************************************/
4195 // Auxiliary functions
4196 Double_t CbmTofDigitize::ComputeClusterAreaOnChannel(Int_t iChanId,
4197  Double_t dClustRadius,
4198  Double_t dClustCentX,
4199  Double_t dClustCentY) {
4200  Double_t dCornersX[4]; // Corners X position
4201  Double_t dCornersY[4]; // Corners X position
4202  Double_t dCornersR[4]; // Distance from cluster center to the corners
4203 
4204  fChannelInfo = fDigiPar->GetCell(iChanId);
4205  Double_t dChanCentPosX = 0.; //fChannelInfo->GetX();
4206  Double_t dChanCentPosY = 0.; //fChannelInfo->GetY();
4207  Double_t dEdgeCentDistX = fChannelInfo->GetSizex() / 2.0;
4208  Double_t dEdgeCentDistY = fChannelInfo->GetSizey() / 2.0;
4209 
4210  // Upper Left
4211  dCornersX[0] = dChanCentPosX - dEdgeCentDistX;
4212  dCornersY[0] = dChanCentPosY + dEdgeCentDistY;
4213  dCornersR[0] = TMath::Sqrt(TMath::Power(dCornersX[0] - dClustCentX, 2)
4214  + TMath::Power(dCornersY[0] - dClustCentY, 2));
4215  // Upper Right
4216  dCornersX[1] = dChanCentPosX + dEdgeCentDistX;
4217  dCornersY[1] = dChanCentPosY + dEdgeCentDistY;
4218  dCornersR[1] = TMath::Sqrt(TMath::Power(dCornersX[1] - dClustCentX, 2)
4219  + TMath::Power(dCornersY[1] - dClustCentY, 2));
4220  // Bottom Right
4221  dCornersX[2] = dChanCentPosX + dEdgeCentDistX;
4222  dCornersY[2] = dChanCentPosY - dEdgeCentDistY;
4223  dCornersR[2] = TMath::Sqrt(TMath::Power(dCornersX[2] - dClustCentX, 2)
4224  + TMath::Power(dCornersY[2] - dClustCentY, 2));
4225  // Bottom Left
4226  dCornersX[3] = dChanCentPosX - dEdgeCentDistX;
4227  dCornersY[3] = dChanCentPosY - dEdgeCentDistY;
4228  dCornersR[3] = TMath::Sqrt(TMath::Power(dCornersX[3] - dClustCentX, 2)
4229  + TMath::Power(dCornersY[3] - dClustCentY, 2));
4230 
4231  Int_t iNbCornersIn = (dCornersR[0] < dClustRadius ? 1 : 0)
4232  + (dCornersR[1] < dClustRadius ? 1 : 0)
4233  + (dCornersR[2] < dClustRadius ? 1 : 0)
4234  + (dCornersR[3] < dClustRadius ? 1 : 0);
4235 
4236  LOG(debug3) << "CbmTofDigitize::ComputeClusterAreaOnChannel => Check "
4237  << iNbCornersIn << " Radius " << dClustRadius << " "
4238  << dCornersR[0] << " " << dCornersR[1] << " " << dCornersR[2]
4239  << " " << dCornersR[3];
4240 
4241  switch (iNbCornersIn) {
4242  case 0: {
4243  // No corner inside the circle
4244  // => if cluster center inside channel: 0-4 disc sections sticking out
4245  // => if cluster center outside channel: 0-1 disc section sticking in
4246  Double_t dEdgeR[4];
4247  dEdgeR[0] =
4248  dChanCentPosX - dEdgeCentDistX - dClustCentX; // Cluster -> Left edge
4249  dEdgeR[1] =
4250  dChanCentPosX + dEdgeCentDistX - dClustCentX; // Cluster -> Right edge
4251  dEdgeR[2] =
4252  dChanCentPosY - dEdgeCentDistY - dClustCentY; // Cluster -> Lower edge
4253  dEdgeR[3] =
4254  dChanCentPosY + dEdgeCentDistY - dClustCentY; // Cluster -> Upper edge
4255  if ((0.0 >= dEdgeR[0]) && (0.0 <= dEdgeR[1]) && (0.0 >= dEdgeR[2])
4256  && (0.0 <= dEdgeR[3])) {
4257  // Cluster center inside channel
4258 
4259  // First disc area
4260  Double_t dArea = dClustRadius * dClustRadius * TMath::Pi();
4261 
4262  LOG(debug3)
4263  << "CbmTofDigitize::ComputeClusterAreaOnChannel => CC in Ch "
4264  << dArea;
4265 
4266  // Now check for each edge if it cuts the cluster circle
4267  // and remove the corresponding disc section if it does
4268  // (no overlap as no corner inside cluster)
4269  if (TMath::Abs(dEdgeR[0]) < dClustRadius)
4270  dArea -= DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR[0]));
4271  if (TMath::Abs(dEdgeR[1]) < dClustRadius)
4272  dArea -= DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR[1]));
4273  if (TMath::Abs(dEdgeR[2]) < dClustRadius)
4274  dArea -= DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR[2]));
4275  if (TMath::Abs(dEdgeR[3]) < dClustRadius)
4276  dArea -= DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR[3]));
4277  if (dArea < 0.0)
4278  LOG(error)
4279  << "CbmTofDigitize::ComputeClusterAreaOnChannel => Neg. area! "
4280  << dArea << " Radius " << dClustRadius << " " << dEdgeR[0] << " "
4281  << dEdgeR[1] << " " << dEdgeR[2] << " " << dEdgeR[3];
4282  return dArea;
4283  } // Cluster center inside channel
4284  else {
4285  // Cluster center outside channel
4286  // At max. one of the edges can cur the cluster circle
4287  // If it is the case, the area on the channel is the disc section area
4288  // If no crossing => no common area of cluster and channel
4289 
4290  LOG(debug3)
4291  << "CbmTofDigitize::ComputeClusterAreaOnChannel => CC out Ch "
4292  << dClustRadius << " " << dEdgeR[0] << " " << dEdgeR[1] << " "
4293  << dEdgeR[2] << " " << dEdgeR[3];
4294 
4295  if (TMath::Abs(dEdgeR[0]) < dClustRadius)
4296  return DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR[0]));
4297  else if (TMath::Abs(dEdgeR[1]) < dClustRadius)
4298  return DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR[1]));
4299  else if (TMath::Abs(dEdgeR[2]) < dClustRadius)
4300  return DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR[2]));
4301  else if (TMath::Abs(dEdgeR[3]) < dClustRadius)
4302  return DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR[3]));
4303  else
4304  return 0.0;
4305  } // Cluster center outside channel
4306  break;
4307  } // case 0
4308  case 1: {
4309  // 1 corner inside the circle
4310  // => we get a "slice" of the cluster disc
4311  // There are then 2 intersection points with the channel, one on each edge
4312  // starting at the included corner. Those 2 points and the corner form a
4313  // triangle. The cluster/channel intersection area is this triangle area +
4314  // the area of the disc section having the 2 intersections for base
4315  Double_t dIntPtX[2] = {0., 0.};
4316  Double_t dIntPtY[2] = {0., 0.};
4317  Double_t dArea = 0.;
4318 
4319  if (dCornersR[0] < dClustRadius) {
4320  // Upper Left corner inside
4321  // Intersection point on left edge
4322  dIntPtX[0] = dCornersX[0];
4323  dIntPtY[0] = CircleIntersectPosY(
4324  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4325  // Intersection point on upper edge
4326  dIntPtX[1] = CircleIntersectPosX(
4327  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4328  dIntPtY[1] = dCornersY[0];
4329 
4330  // First triangle area
4331  dArea = TriangleArea(dCornersX[0],
4332  dCornersY[0],
4333  dIntPtX[0],
4334  dIntPtY[0],
4335  dIntPtX[1],
4336  dIntPtY[1]);
4337  } else if (dCornersR[1] < dClustRadius) {
4338  // Upper Right corner inside
4339  // Intersection point on Right edge
4340  dIntPtX[0] = dCornersX[1];
4341  dIntPtY[0] = CircleIntersectPosY(
4342  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4343  // Intersection point on upper edge
4344  dIntPtX[1] = CircleIntersectPosX(
4345  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4346  dIntPtY[1] = dCornersY[1];
4347 
4348  // First triangle area
4349  dArea = TriangleArea(dCornersX[1],
4350  dCornersY[1],
4351  dIntPtX[0],
4352  dIntPtY[0],
4353  dIntPtX[1],
4354  dIntPtY[1]);
4355  } else if (dCornersR[2] < dClustRadius) {
4356  // Bottom Right corner inside
4357  // Intersection point on Right edge
4358  dIntPtX[0] = dCornersX[2];
4359  dIntPtY[0] = CircleIntersectPosY(
4360  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4361  // Intersection point on lower edge
4362  dIntPtX[1] = CircleIntersectPosX(
4363  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4364  dIntPtY[1] = dCornersY[2];
4365 
4366  // First triangle area
4367  dArea = TriangleArea(dCornersX[2],
4368  dCornersY[2],
4369  dIntPtX[0],
4370  dIntPtY[0],
4371  dIntPtX[1],
4372  dIntPtY[1]);
4373  } else if (dCornersR[3] < dClustRadius) {
4374  // Bottom Left corner inside
4375  // Intersection point on left edge
4376  dIntPtX[0] = dCornersX[3];
4377  dIntPtY[0] = CircleIntersectPosY(
4378  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4379  // Intersection point on lower edge
4380  dIntPtX[1] = CircleIntersectPosX(
4381  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4382  dIntPtY[1] = dCornersY[3];
4383 
4384  // First triangle area
4385  dArea = TriangleArea(dCornersX[3],
4386  dCornersY[3],
4387  dIntPtX[0],
4388  dIntPtY[0],
4389  dIntPtX[1],
4390  dIntPtY[1]);
4391  }
4392 
4393  // Then add disc section area
4394  // Compute the distance between base and cluster center
4395  Double_t dSecBaseD = DistanceCircleToBase(
4396  dClustRadius, dIntPtX[0], dIntPtY[0], dIntPtX[1], dIntPtY[1]);
4397  // Add computed area
4398  dArea += DiscSectionArea(dClustRadius, dSecBaseD);
4399 
4400  return dArea;
4401  break;
4402  } // case 1
4403  case 2: {
4404  // 2 corners inside the circle
4405  // => 1 edge is fully included
4406  // Each of the edges in the other direction has 1 intersection point
4407  // with the circle. The area between the included edge and the line
4408  // they form can be otained by summing the area of 2 triangles, using
4409  // 1 of the corners inside the circle.
4410  // For the disc section having these points for base, there are 2
4411  // cases. Either it is fully contained inside the channel, or it sticks
4412  // out, making 2 intersection points on the edge opposite to the
4413  // included one. In the second case the area of this second disc section
4414  // has to be subtracted.
4415  Bool_t bFirstCorner = kTRUE;
4416  Int_t iCorners[2];
4417 
4418  iCorners[0] = -1;
4419  iCorners[1] = -1;
4420  if (dCornersR[0] < dClustRadius) {
4421  // Upper Left corner inside
4422  iCorners[0] = 0;
4423  bFirstCorner = kFALSE;
4424  } // if( dCornersR[0] < dClustRadius )
4425  if (dCornersR[1] < dClustRadius) {
4426  // Upper Right corner inside
4427  if (kTRUE == bFirstCorner) {
4428  iCorners[0] = 1;
4429  bFirstCorner = kFALSE;
4430  } // if( kTRUE == bFirstCorner )
4431  else
4432  iCorners[1] = 1;
4433  } // else if( dCornersR[1] < dClustRadius )
4434  if (dCornersR[2] < dClustRadius) {
4435  // Bottom Right corner inside
4436  if (kTRUE == bFirstCorner) {
4437  iCorners[0] = 2;
4438  bFirstCorner = kFALSE;
4439  } // if( kTRUE == bFirstCorner )
4440  else
4441  iCorners[1] = 2;
4442  } // else if( dCornersR[2] < dClustRadius )
4443  if (dCornersR[3] < dClustRadius) {
4444  // Bottom Left corner inside
4445  // Has to be the 2nd one if there
4446  iCorners[1] = 3;
4447  } // else if( dCornersR[3] < dClustRadius )
4448 
4449  LOG(debug3) << "Corners In check " << (iCorners[0]) << " "
4450  << (iCorners[1]);
4451  // Get the interesting points info
4452  Double_t dIntPtX[2] = {0., 0.};
4453  Double_t dIntPtY[2] = {0., 0.};
4454  Double_t dEdgeR = 0.;
4455  Int_t iOppCorner = 0;
4456  if (0 == iCorners[0] && 1 == iCorners[1]) {
4457  LOG(debug3) << "Test upper edge";
4458  // Upper edge
4459  // 1st Intersection point on left edge
4460  dIntPtX[0] = dCornersX[0];
4461  dIntPtY[0] = CircleIntersectPosY(
4462  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4463  // 2nd Intersection point on right edge
4464  dIntPtX[1] = dCornersX[1];
4465  dIntPtY[1] = CircleIntersectPosY(
4466  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4467  // Corner opposed to first intersection
4468  iOppCorner = 1; // <- Upper Right corner
4469  // Distance between out edge and cluster center
4470  dEdgeR = dChanCentPosY - dEdgeCentDistY
4471  - dClustCentY; // Cluster -> Lower edge
4472  } else if (1 == iCorners[0] && 2 == iCorners[1]) {
4473  LOG(debug3) << "Test right edge";
4474  // Right edge
4475  // Intersection point on upper edge
4476  dIntPtX[0] = CircleIntersectPosX(
4477  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4478  dIntPtY[0] = dCornersY[1];
4479  // Intersection point on lower edge
4480  dIntPtX[1] = CircleIntersectPosX(
4481  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4482  dIntPtY[1] = dCornersY[2];
4483  // Corner opposed to first intersection
4484  iOppCorner = 2; // <- Bottom Right corner
4485  // Distance between out edge and cluster center
4486  dEdgeR =
4487  dChanCentPosX - dEdgeCentDistX - dClustCentX; // Cluster -> Left edge
4488  } else if (2 == iCorners[0] && 3 == iCorners[1]) {
4489  LOG(debug3) << "Test lower edge";
4490  // Lower edge
4491  // 1st Intersection point on right edge
4492  dIntPtX[0] = dCornersX[2];
4493  dIntPtY[0] = CircleIntersectPosY(
4494  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4495  // 2nd Intersection point on left edge
4496  dIntPtX[1] = dCornersX[3];
4497  dIntPtY[1] = CircleIntersectPosY(
4498  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4499  // Corner opposed to first intersection
4500  iOppCorner = 3; // <- Bottom left corner
4501  // Distance between out edge and cluster center
4502  dEdgeR = dChanCentPosY + dEdgeCentDistY
4503  - dClustCentY; // Cluster -> Upper edge
4504  } else if (0 == iCorners[0] && 3 == iCorners[1]) {
4505  LOG(debug3) << "Test left edge";
4506  // Left edge
4507  // Intersection point on lower edge
4508  dIntPtX[0] = CircleIntersectPosX(
4509  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4510  dIntPtY[0] = dCornersY[3];
4511  // Intersection point on upper edge
4512  dIntPtX[1] = CircleIntersectPosX(
4513  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4514  dIntPtY[1] = dCornersY[0];
4515  // Corner opposed to first intersection
4516  iOppCorner = 0; // <- Upper left corner
4517  // Distance between out edge and cluster center
4518  dEdgeR = dChanCentPosX + dEdgeCentDistX
4519  - dClustCentX; // Cluster -> Right edge
4520  }
4521 
4522  // First triangle: The 2 corners and the 1st intersection
4523  Double_t dArea = TriangleArea(dCornersX[iCorners[0]],
4524  dCornersY[iCorners[0]],
4525  dCornersX[iCorners[1]],
4526  dCornersY[iCorners[1]],
4527  dIntPtX[0],
4528  dIntPtY[0]);
4529  // Second triangle: The corners opposed to first intersection and the 2 intersections
4530  dArea += TriangleArea(dCornersX[iOppCorner],
4531  dCornersY[iOppCorner],
4532  dIntPtX[0],
4533  dIntPtY[0],
4534  dIntPtX[1],
4535  dIntPtY[1]);
4536 
4537  // Now add the disc section
4538  // Compute the distance between base and cluster center
4539  Double_t dSecBaseD = DistanceCircleToBase(
4540  dClustRadius, dIntPtX[0], dIntPtY[0], dIntPtX[1], dIntPtY[1]);
4541  // Add computed area
4542  dArea += DiscSectionArea(dClustRadius, dSecBaseD);
4543 
4544  // Use the distance between the cluster center and the opposite edge to
4545  // check if we need to remove a part sticking out
4546 
4547  // Now check for this edge if it cuts the cluster circle
4548  // and remove the corresponding disc section if it does
4549  if (TMath::Abs(dEdgeR) < dClustRadius)
4550  dArea -= DiscSectionArea(dClustRadius, TMath::Abs(dEdgeR));
4551 
4552  return dArea;
4553  break;
4554  } // case 2
4555  case 3: {
4556  // 3 corners inside the circle
4557  // => 2 edges fully included, 1 intersection on each of the 2 others
4558  // In this case the ovelapped area is equal to the full channel area
4559  // minus the area of the triangle formed by the out corner and the 2
4560  // intersection, plus the area of the disc section having the 2
4561  // intersection points for base.
4562  Int_t iOutCorner = 0;
4563  Double_t dIntPtX[2] = {0., 0.};
4564  Double_t dIntPtY[2] = {0., 0.};
4565 
4566  if (dCornersR[0] > dClustRadius) {
4567  // Upper Left corner out
4568  iOutCorner = 0;
4569  // Intersection on upper edge
4570  dIntPtX[0] = CircleIntersectPosX(
4571  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4572  dIntPtY[0] = dCornersY[0];
4573  // Intersection on left edge
4574  dIntPtX[1] = dCornersX[0];
4575  dIntPtY[1] = CircleIntersectPosY(
4576  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4577  } // if( dCornersR[0] > dClustRadius )
4578  else if (dCornersR[1] > dClustRadius) {
4579  // Upper Right corner out
4580  iOutCorner = 1;
4581  // Intersection on upper edge
4582  dIntPtX[0] = CircleIntersectPosX(
4583  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4584  dIntPtY[0] = dCornersY[1];
4585  // Intersection on right edge
4586  dIntPtX[1] = dCornersX[1];
4587  dIntPtY[1] = CircleIntersectPosY(
4588  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4589  } // else if( dCornersR[1] > dClustRadius )
4590  else if (dCornersR[2] > dClustRadius) {
4591  // Bottom Right corner out
4592  iOutCorner = 2;
4593  // Intersection on bottom edge
4594  dIntPtX[0] = CircleIntersectPosX(
4595  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4596  dIntPtY[0] = dCornersY[2];
4597  // Intersection on right edge
4598  dIntPtX[1] = dCornersX[2];
4599  dIntPtY[1] = CircleIntersectPosY(
4600  iChanId, dClustRadius, dClustCentX, dClustCentY, kTRUE);
4601  } // else if( dCornersR[2] > dClustRadius )
4602  else if (dCornersR[3] > dClustRadius) {
4603  // Bottom Left corner out
4604  iOutCorner = 3;
4605  // Intersection on bottom edge
4606  dIntPtX[0] = CircleIntersectPosX(
4607  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4608  dIntPtY[0] = dCornersY[3];
4609  // Intersection on left edge
4610  dIntPtX[1] = dCornersX[3];
4611  dIntPtY[1] = CircleIntersectPosY(
4612  iChanId, dClustRadius, dClustCentX, dClustCentY, kFALSE);
4613  } // else if( dCornersR[3] > dClustRadius )
4614 
4615  // First get full channel area
4616  Double_t dArea = (fChannelInfo->GetSizex()) * (fChannelInfo->GetSizey());
4617 
4618  // Then subtract the triangle area
4619  dArea -= TriangleArea(dCornersX[iOutCorner],
4620  dCornersY[iOutCorner],
4621  dIntPtX[0],
4622  dIntPtY[0],
4623  dIntPtX[1],
4624  dIntPtY[1]);
4625 
4626  // Finally add the disc section area
4627  // Compute the distance between base and cluster center
4628  Double_t dSecBaseD = DistanceCircleToBase(
4629  dClustRadius, dIntPtX[0], dIntPtY[0], dIntPtX[1], dIntPtY[1]);
4630  // Add computed area
4631  dArea += DiscSectionArea(dClustRadius, dSecBaseD);
4632 
4633  return dArea;
4634  break;
4635  } // case 3
4636  case 4: {
4637  // All in circle => channel fully contained!
4638  // => Area of Cluster/channel intersection = channel area
4639  return (fChannelInfo->GetSizex()) * (fChannelInfo->GetSizey());
4640  break;
4641  } // case 4
4642  default:
4643  // Should never happen !!!
4644  return 0;
4645  break;
4646  } // switch( iNbCornersIn )
4647 }
4648 /************************************************************************************/
4649 // Auxiliary functions
4650 // Area
4651 Double_t CbmTofDigitize::TriangleArea(Double_t dXa,
4652  Double_t dYa,
4653  Double_t dXb,
4654  Double_t dYb,
4655  Double_t dXc,
4656  Double_t dYc) {
4657  Double_t dArea =
4658  0.5 * ((dXa - dXc) * (dYb - dYa) - (dXa - dXb) * (dYc - dYa));
4659  return TMath::Abs(dArea);
4660 }
4661 Double_t CbmTofDigitize::DiscSectionArea(Double_t dDiscRadius,
4662  Double_t dDistBaseToCenter) {
4663  if (dDiscRadius < dDistBaseToCenter || dDiscRadius <= 0
4664  || dDistBaseToCenter < 0) {
4665  LOG(error) << "CbmTofDigitize::DiscSectionArea => Invalid Input values!! "
4666  "Disc radius "
4667  << dDiscRadius << " and Base distance " << dDistBaseToCenter;
4668  return 0.0;
4669  }
4670  // First Half disc area
4671  Double_t dArea = dDiscRadius * dDiscRadius * TMath::Pi() / 2.0;
4672 
4673  // Then remove the "rectangle" middle part of the area between base and center
4674  dArea -= dDistBaseToCenter
4675  * TMath::Sqrt(dDiscRadius * dDiscRadius
4676  - dDistBaseToCenter * dDistBaseToCenter);
4677  // Finally remove the "arc like" side parts of the area between base and center
4678  dArea -=
4679  dDiscRadius * dDiscRadius * TMath::ASin(dDistBaseToCenter / dDiscRadius);
4680 
4681  return dArea;
4682 }
4683 //----------------------------------------------------------------------------//
4684 // Points
4685 Double_t CbmTofDigitize::CircleIntersectPosX(Int_t iChanId,
4686  Double_t dClustRadius,
4687  Double_t dClustCentX,
4688  Double_t dClustCentY,
4689  Bool_t bUpperSide) {
4690  fChannelInfo = fDigiPar->GetCell(iChanId);
4691  Double_t dChanCentPosX = 0.; //fChannelInfo->GetX();
4692  Double_t dChanCentPosY = 0.; //fChannelInfo->GetY();
4693  Double_t dEdgeCentDistX = fChannelInfo->GetSizex() / 2.0;
4694  Double_t dEdgeCentDistY = fChannelInfo->GetSizey() / 2.0;
4695 
4696  Double_t dEdgePosY =
4697  dChanCentPosY + (kTRUE == bUpperSide ? 1 : -1) * dEdgeCentDistY;
4698 
4699  if (dChanCentPosX == dClustCentX) {
4700  // Clustered centered on channel center
4701  // => a single intersection means that the distance between cluster and edge is exactly
4702  // the radius
4703  if (TMath::Abs(dEdgePosY - dClustCentY) == dClustRadius)
4704  // => Intersection at same X position as channel center
4705  return dChanCentPosX;
4706  // Other values mean 0 or 2 intersections in edge range
4707  // => return 0.0, faulty call to this function
4708  else {
4709  LOG(error) << "CbmTofDigitize::CircleIntersectPosX => Invalid values: "
4710  << " 0 or 2 intersections with cluster aligned on channel ";
4711  return 0.0;
4712  } // else of if( dEdgeCentDistY == dClustRadius )
4713  } // if( dChanCentPosX == dClustCentX )
4714  else {
4715  Double_t dRoot =
4716  dClustRadius * dClustRadius - TMath::Power(dClustCentY - dEdgePosY, 2);
4717  if (0.0 <= dRoot) {
4718  // Null root means 1 single solution
4719  // Positive root means 2 solutions
4720  /*
4721  // If the circle center is more to the left than the channel center and
4722  // there are 2 solutions, assuming only 1 intersection exists inside the
4723  // edge boundaries, it will be the rightmost one.
4724  // If the circle center is more to the left than the channel center, it
4725  // will be the leftmost one.
4726  */
4727  Double_t dDistX = TMath::Sqrt(dRoot);
4728  Double_t dSign =
4729  ((dClustCentX - dDistX > dChanCentPosX - dEdgeCentDistX)
4730  && (dClustCentX - dDistX < dChanCentPosX + dEdgeCentDistX)
4731  ? -1.0
4732  : 1.0);
4733  return dClustCentX + dSign * dDistX;
4734  } // if( 0.0 <= dRoot )
4735  // Negative root means no intersection at all between
4736  // the circle and the line generated by this edge!!
4737  // => return 0.0, faulty call to this function
4738  else {
4739  LOG(error) << "CbmTofDigitize::CircleIntersectPosX => Invalid values: "
4740  << " 0 intersection at all (negative root = " << dRoot << ") ";
4741  TString sTemp =
4742  Form(" Radius = %5.4f ClustX = %6.3f ClustY = %6.3f Side = %1d EdgeY = "
4743  "%6.3f ChanX = %6.3f ChanY = %6.3f Channel = %6d",
4744  dClustRadius,
4745  dClustCentX,
4746  dClustCentY,
4747  bUpperSide,
4748  dEdgePosY,
4749  dChanCentPosX,
4750  dChanCentPosY,
4751  iChanId);
4752  LOG(error) << "CbmTofDigitize::CircleIntersectPosX => Input values: "
4753  << sTemp;
4754  return 0.0;
4755  } // else of if( 0.0 <= dRoot )
4756  } // else of if( dChanCentPosX == dClustCentX )
4757 }
4758 Double_t CbmTofDigitize::CircleIntersectPosY(Int_t iChanId,
4759  Double_t dClustRadius,
4760  Double_t dClustCentX,
4761  Double_t dClustCentY,
4762  Bool_t bRightSide) {
4763  fChannelInfo = fDigiPar->GetCell(iChanId);
4764  Double_t dChanCentPosX = 0.; //fChannelInfo->GetX();
4765  Double_t dChanCentPosY = 0.; //fChannelInfo->GetY();
4766  Double_t dEdgeCentDistX = fChannelInfo->GetSizex() / 2.0;
4767  Double_t dEdgeCentDistY = fChannelInfo->GetSizey() / 2.0;
4768 
4769  Double_t dEdgePosX =
4770  dChanCentPosX + (kTRUE == bRightSide ? 1 : -1) * dEdgeCentDistX;
4771 
4772  if (dChanCentPosY == dClustCentY) {
4773  // Clustered centered on channel center
4774  // => a single intersection means that the distance between cluster and edge is exactly
4775  // the radius
4776  if (TMath::Abs(dEdgePosX - dClustCentX) == dClustRadius)
4777  // => Intersection at same X position as channel center
4778  return dChanCentPosY;
4779  // Other values mean 0 or 2 intersections in edge range
4780  // => return 0.0, faulty call to this function
4781  else {
4782  LOG(error) << "CbmTofDigitize::CircleIntersectPosY => Invalid values: "
4783  << " 0 or 2 intersections with cluster aligned on channel ";
4784  return 0.0;
4785  } // else of if( dEdgeCentDistX == dClustRadius )
4786  } // if( dChanCentPosX == dClustCentX )
4787  else {
4788  Double_t dRoot =
4789  dClustRadius * dClustRadius - TMath::Power(dClustCentX - dEdgePosX, 2);
4790  if (0.0 <= dRoot) {
4791  // Null root means 1 single solution
4792  // Positive root means 2 solutions
4793  /*
4794  // If the circle center is higher than the channel center and
4795  // there are 2 solutions, assuming only 1 intersection exists inside the
4796  // edge boundaries, it will be lower one.
4797  // If the circle center is more to the left than the channel center, it
4798  // will be the higher one.
4799  Double_t dSign = ( dClustCentY > dChanCentPosY ? -1.0 : 1.0);
4800 
4801  return dChanCentPosY + dSign * TMath::Sqrt( dRoot );
4802  */
4803  Double_t dDistY = TMath::Sqrt(dRoot);
4804  Double_t dSign =
4805  ((dClustCentY - dDistY > dChanCentPosY - dEdgeCentDistY)
4806  && (dClustCentY - dDistY < dChanCentPosY + dEdgeCentDistY)
4807  ? -1.0
4808  : 1.0);
4809  return dClustCentY + dSign * dDistY;
4810  } // if( 0.0 <= dRoot )
4811  // Negative root means no intersection at all between
4812  // the circle and the line generated by this edge!!
4813  // => return 0.0, faulty call to this function
4814  else {
4815  LOG(error) << "CbmTofDigitize::CircleIntersectPosY => Invalid values: "
4816  << " 0 intersection at all (negative root = " << dRoot << ") ";
4817  TString sTemp =
4818  Form(" Radius = %5.4f ClustX = %6.3f ClustY = %6.3f Side = %1d EdgeX = "
4819  "%6.3f ChanX = %6.3f ChanY = %6.3f Channel = %6d",
4820  dClustRadius,
4821  dClustCentX,
4822  dClustCentY,
4823  bRightSide,
4824  dEdgePosX,
4825  dChanCentPosX,
4826  dChanCentPosY,
4827  iChanId);
4828  LOG(error) << "CbmTofDigitize::CircleIntersectPosY => Input values: "
4829  << sTemp;
4830  return 0.0;
4831  } // else of if( 0.0 <= dRoot )
4832  } // else of if( dChanCentPosX == dClustCentX )
4833 }
4834 Double_t CbmTofDigitize::DistanceCircleToBase(Double_t dClustRadius,
4835  Double_t dBaseXa,
4836  Double_t dBaseYa,
4837  Double_t dBaseXb,
4838  Double_t dBaseYb) {
4839  // Both base endpoints should be on the circle, forming an isoscele triangle
4840  // with the circle center.
4841  // The distance to the base is then the height of the isoscele triangle.
4842  Double_t dBaseLength = TMath::Sqrt(TMath::Power(dBaseXb - dBaseXa, 2)
4843  + TMath::Power(dBaseYb - dBaseYa, 2));
4844  Double_t dRoot = dClustRadius * dClustRadius - dBaseLength * dBaseLength / 4;
4845  if (0.0 <= dRoot)
4846  return TMath::Sqrt(dRoot);
4847  else {
4848  LOG(error) << "CbmTofDigitize::DistanceCircleToBase => Invalid values: "
4849  << " base end-points not on circle (negative root" << dRoot
4850  << ") ";
4851  TString sTemp =
4852  Form(" Radius = %5.4f Xa = %6.3f Ya = %6.3f Xb = %6.3f Yb = %6.3f ",
4853  dClustRadius,
4854  dBaseXa,
4855  dBaseYa,
4856  dBaseXb,
4857  dBaseYb);
4858  LOG(error) << "CbmTofDigitize::DistanceCircleToBase => Input values: "
4859  << sTemp;
4860  return 0.0;
4861  } // else of if( 0.0 <= dRoot )
4862 }
4863 
4864 Bool_t CbmTofDigitize::CompareTimes(CbmTofDigi* pDigi1, CbmTofDigi* pDigi2) {
4865  return (pDigi1->GetTime() < pDigi2->GetTime());
4866 }
4867 
4868 
4895 /************************************************************************************/
CbmTofDigitize::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: CbmTofDigitize.cxx:271
CbmTofDigiBdfPar.h
CbmTofCell::GetZ
Double_t GetZ() const
Definition: CbmTofCell.h:38
CbmTofDigitize::TriangleArea
Double_t TriangleArea(Double_t dXa, Double_t dYa, Double_t dXb, Double_t dYb, Double_t dXc, Double_t dYc)
Compute triangle area from its corners.
Definition: CbmTofDigitize.cxx:4651
CbmTofDigitize::FillHistos
Bool_t FillHistos()
Definition: CbmTofDigitize.cxx:880
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmTofDigitize.h
CbmTofDetectorId::GetGap
virtual Int_t GetGap(const Int_t detectorId)=0
CbmTofGeoHandler::GetGeoVersion
Int_t GetGeoVersion()
Definition: CbmTofGeoHandler.h:45
CbmTofDigitize::fStop
TTimeStamp fStop
Definition: CbmTofDigitize.h:332
CbmMatch
Definition: CbmMatch.h:22
CbmTofDigiBdfPar::printParams
void printParams()
Definition: CbmTofDigiBdfPar.cxx:667
CbmTofDigitize::fDigiBdfPar
CbmTofDigiBdfPar * fDigiBdfPar
Definition: CbmTofDigitize.h:267
CbmTofCell::GetSizex
Double_t GetSizex() const
Definition: CbmTofCell.h:40
CbmTofDigitize::fbTimeBasedOutput
Bool_t fbTimeBasedOutput
Definition: CbmTofDigitize.h:347
CbmTofDigitize::fdNofPointsTot
Double_t fdNofPointsTot
Total number of points processed.
Definition: CbmTofDigitize.h:339
CbmTofDigitize::fhDigiTime
TH1 * fhDigiTime
Definition: CbmTofDigitize.h:319
CbmTofDigitize::fdDigiTimeConvFactor
Double_t fdDigiTimeConvFactor
Definition: CbmTofDigitize.h:353
CbmTofGeoHandler::GetCounter
Int_t GetCounter(Int_t uniqueId)
Definition: CbmTofGeoHandler.cxx:469
CbmMatch::GetLink
const CbmLink & GetLink(Int_t i) const
Definition: CbmMatch.h:35
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
CbmTofDigitize::MergeSameChanDigis
Bool_t MergeSameChanDigis()
Merge the digis on he same readout channel.
Definition: CbmTofDigitize.cxx:1143
CbmTofDigitize::RegisterInputs
Bool_t RegisterInputs()
Recover pointer on input TClonesArray: TofPoints, ...
Definition: CbmTofDigitize.cxx:351
CbmTofDigiBdfPar::GetClustSizeHist
TH1 * GetClustSizeHist(Int_t iSmType) const
Definition: CbmTofDigiBdfPar.cxx:637
CbmTofDigi::SetTime
void SetTime(Double_t time)
Definition: CbmTofDigi.h:153
CbmDigitizeBase::GetEventInfo
void GetEventInfo()
Get event information.
Definition: CbmDigitizeBase.cxx:48
CbmTofDigitize::fhTofPtsInTrkVsGapInd
TH2 * fhTofPtsInTrkVsGapInd
Definition: CbmTofDigitize.h:301
CbmTofDigiBdfPar::GetResolution
Double_t GetResolution(Int_t iSmType) const
Definition: CbmTofDigiBdfPar.cxx:623
CbmTofDigiBdfPar::UseOneGapPerTrk
Bool_t UseOneGapPerTrk() const
Definition: CbmTofDigiBdfPar.h:89
CompTimesExp
Definition: CbmTofDigitize.cxx:63
CbmTofDigitize::fhTofPtsInTrkVsGapIndPrm
TH2 * fhTofPtsInTrkVsGapIndPrm
Definition: CbmTofDigitize.h:302
CbmTofDigitize::fdDigitizeTime
Double_t fdDigitizeTime
Definition: CbmTofDigitize.h:333
CbmTofDetectorInfo::fCell
Int_t fCell
Definition: CbmTofDetectorId.h:62
CbmTofDetectorId::GetSMType
virtual Int_t GetSMType(const Int_t detectorId)=0
CbmTofDigitize::GenerateClusterRadius
Double_t GenerateClusterRadius(Int_t iSmType, Int_t iRpc)
Generate a value for the cluster radius from the beamtime data corresponding.
Definition: CbmTofDigitize.cxx:1409
CbmTofDigiBdfPar::GetNbChan
Int_t GetNbChan(Int_t iSmType, Int_t iRpc) const
Definition: CbmTofDigiBdfPar.cxx:570
CbmTofDigiPar::GetNrOfModules
Int_t GetNrOfModules()
Definition: CbmTofDigiPar.h:44
CbmTofDigitize::fvdSignalVelocityRpc
std::vector< std::vector< std::vector< Double_t > > > fvdSignalVelocityRpc
Definition: CbmTofDigitize.h:256
CbmTofDigitize::DeleteHistos
Bool_t DeleteHistos()
Definition: CbmTofDigitize.cxx:1096
CbmTofDigiBdfPar::GetGapEfficiency
Double_t GetGapEfficiency(Int_t iSmType, Int_t iRpc) const
Definition: CbmTofDigiBdfPar.cxx:613
CbmTofDigitize::fhMeanDigiPerTrack
TH1 * fhMeanDigiPerTrack
Definition: CbmTofDigitize.h:316
CbmTofDigitize::fh1ClusterSizeProb
std::vector< TH1 * > fh1ClusterSizeProb
Definition: CbmTofDigitize.h:252
CbmTofDigitize::SetInputFileName
void SetInputFileName(TString FileName)
Definition: CbmTofDigitize.h:73
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmTofDigitize::fdTimeTot
Double_t fdTimeTot
Total execution time.
Definition: CbmTofDigitize.h:341
CbmTofDigiBdfPar::GetNbRpc
Int_t GetNbRpc(Int_t iSmType) const
Definition: CbmTofDigiBdfPar.cxx:519
CbmTofDigitize::fhElecChOccup
TH1 * fhElecChOccup
Definition: CbmTofDigitize.h:324
CbmTofDigitize::fdFeeTotThr
Double_t fdFeeTotThr
Definition: CbmTofDigitize.h:247
CbmTofDetectorId::GetCounter
virtual Int_t GetCounter(const Int_t detectorId)=0
CbmTofDigitize::fvRpcChOffs
std::vector< std::vector< std::vector< Int_t > > > fvRpcChOffs
Definition: CbmTofDigitize.h:270
CbmTofDigitize::~CbmTofDigitize
virtual ~CbmTofDigitize()
Destructor.
Definition: CbmTofDigitize.cxx:223
CbmTofDigitize::fhEvtProcTime
TH1 * fhEvtProcTime
Definition: CbmTofDigitize.h:312
CbmDigitizeBase::fCurrentMCEntry
Int_t fCurrentMCEntry
Number of current MC event.
Definition: CbmDigitizeBase.h:164
CbmTofGeoHandler::GetGap
Int_t GetGap(Int_t uniqueId)
Definition: CbmTofGeoHandler.cxx:474
CbmTofDetectorId_v12b
Definition: CbmTofDetectorId_v12b.h:33
CbmTofDigiBdfPar::GetChanOrient
Int_t GetChanOrient(Int_t iSmType, Int_t iRpc) const
Definition: CbmTofDigiBdfPar.cxx:590
CbmTofDigiBdfPar::GetSigVel
Double_t GetSigVel(Int_t iSmType, Int_t iSm, Int_t iRpc) const
Definition: CbmTofDigiBdfPar.cxx:546
CbmTofDigiBdfPar::GetNbSmTypes
Int_t GetNbSmTypes() const
Definition: CbmTofDigiBdfPar.h:56
CbmTofDigitize::fbAllowPointsWithoutTrack
Bool_t fbAllowPointsWithoutTrack
Definition: CbmTofDigitize.h:348
CbmTofDigiBdfPar::GetClusterRadiusModel
Int_t GetClusterRadiusModel() const
Definition: CbmTofDigiBdfPar.h:73
CbmTofDigiBdfPar::GetFeeGainSigma
Double_t GetFeeGainSigma() const
Definition: CbmTofDigiBdfPar.h:47
CbmTofDigitize::fStart
TTimeStamp fStart
Definition: CbmTofDigitize.h:331
CbmTofDigi.h
ECbmModuleId::kTof
@ kTof
Time-of-flight Detector.
CbmTofDigitize::fiNbDigis
Int_t fiNbDigis
Definition: CbmTofDigitize.h:295
CbmTofDigiPar.h
CbmTofDigitize::SetHistoFileName
Bool_t SetHistoFileName(TString sFilenameIn="./tofDigiBdf.hst.root")
Definition: CbmTofDigitize.cxx:1027
CbmMatch.h
CbmTofDigitize::fMcTracksColl
TClonesArray * fMcTracksColl
Definition: CbmTofDigitize.h:274
CbmTofDigiPar::GetCell
CbmTofCell * GetCell(Int_t i)
Definition: CbmTofDigiPar.h:48
CbmTofDigitize::DigitizeGaussCharge
Bool_t DigitizeGaussCharge()
Convert TofPoints in input TClonesArray to Tof Digis using an approximation of the.
Definition: CbmTofDigitize.cxx:3123
CbmTofDigiBdfPar::GetDeadtime
Double_t GetDeadtime() const
Definition: CbmTofDigiBdfPar.h:51
CbmTofGeoHandler::GetSModule
Int_t GetSModule(Int_t uniqueId)
Definition: CbmTofGeoHandler.cxx:464
CbmTofCreateDigiPar.h
CbmTofDigitize::fhPtTime
TH1 * fhPtTime
Definition: CbmTofDigitize.h:318
CbmTofDigitize::fdFeeGainSigma
Double_t fdFeeGainSigma
Definition: CbmTofDigitize.h:246
CbmTofDigitize::SetParContainers
virtual void SetParContainers()
Inherited from FairTask.
Definition: CbmTofDigitize.cxx:259
CbmTofGeoHandler.h
CbmTofDigitize::ComputeClusterAreaOnChannel
Double_t ComputeClusterAreaOnChannel(Int_t iChanId, Double_t dClustRadius, Double_t dClustCentX, Double_t dClustCentY)
Compute geometrical intersection area of a cluster and a channel.
Definition: CbmTofDigitize.cxx:4196
CbmTofDetectorId_v14a.h
CbmTofDigitize::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmTofDigitize.cxx:230
fCurrentMCEntry
Int_t fCurrentMCEntry
Definition: CbmTofDigitize.cxx:82
CbmTofDigitize::fhDigiTimeResB
TH2 * fhDigiTimeResB
Definition: CbmTofDigitize.h:321
TrbNetState::DEBUG
@ DEBUG
CbmTofDigitize::fhNbTracksEvtElCh
TH2 * fhNbTracksEvtElCh
Definition: CbmTofDigitize.h:327
CbmTofDigitize::fbMonitorHistos
Bool_t fbMonitorHistos
Definition: CbmTofDigitize.h:345
CbmTofGeoHandler::GetCellId
Int_t GetCellId(Int_t uniqueId)
Definition: CbmTofGeoHandler.cxx:493
CbmTofDigitize::InitParameters
Bool_t InitParameters()
Initialize other parameters not included in parameter classes.
Definition: CbmTofDigitize.cxx:370
CbmTofDigiBdfPar::GetClusterModel
Int_t GetClusterModel() const
Definition: CbmTofDigiBdfPar.h:92
CbmTofDigiBdfPar::UseOnlyPrimaries
Bool_t UseOnlyPrimaries() const
Definition: CbmTofDigiBdfPar.h:87
CbmTofDigitize::fhProcTimeEvtSize
TH2 * fhProcTimeEvtSize
Definition: CbmTofDigitize.h:315
CbmMCTrack::GetNPoints
Int_t GetNPoints(ECbmModuleId detId) const
Definition: CbmMCTrack.cxx:186
CbmTofGeoHandler
Definition: CbmTofGeoHandler.h:30
CbmTofDigitize::fhTofPtsPosVsGap
TH2 * fhTofPtsPosVsGap[10]
Definition: CbmTofDigitize.h:304
k12b
@ k12b
Definition: CbmTofGeoHandler.h:17
CbmDigitize
Base class template for CBM digitisation tasks.
Definition: CbmDigitize.h:39
CbmTofCell::GetX
Double_t GetX() const
Definition: CbmTofCell.h:36
CbmTofDetectorId_v12b.h
CbmTofDigitize::fhMultiDigiEvtElCh
TH1 * fhMultiDigiEvtElCh
Definition: CbmTofDigitize.h:325
CbmDigitizeBase::fEventMode
Bool_t fEventMode
Definition: CbmDigitizeBase.h:159
CbmTofDigitize::fdNofTofMcTrkTot
Double_t fdNofTofMcTrkTot
Total number of MC tracks with TOF points.
Definition: CbmTofDigitize.h:338
CbmTofDigitize::fGeoHandler
CbmTofGeoHandler * fGeoHandler
Definition: CbmTofDigitize.h:263
CbmTofDigitize::fh1ClusterTotProb
std::vector< TH1 * > fh1ClusterTotProb
Definition: CbmTofDigitize.h:253
CbmTofDigitize::fiNbElecChTot
Int_t fiNbElecChTot
Definition: CbmTofDigitize.h:268
CbmDigitizeBase::fCurrentEvent
Int_t fCurrentEvent
Number of current input.
Definition: CbmDigitizeBase.h:163
k14a
@ k14a
Definition: CbmTofGeoHandler.h:17
CbmMatch::AddLink
void AddLink(const CbmLink &newLink)
Definition: CbmMatch.cxx:42
CbmTofDigitize::fTofPointsColl
TClonesArray * fTofPointsColl
Definition: CbmTofDigitize.h:273
CbmTofDigi::GetTime
Double_t GetTime() const
Inherited from CbmDigi.
Definition: CbmTofDigi.h:111
CbmTofAddress.h
CbmTofDetectorId_v14a
Definition: CbmTofDetectorId_v14a.h:36
CbmTofDigitize::fdMergeTime
Double_t fdMergeTime
Definition: CbmTofDigitize.h:334
CbmTofGeoHandler::Init
Int_t Init(Bool_t isSimulation=kFALSE)
Definition: CbmTofGeoHandler.cxx:39
CbmTofDigitize::fdSignalPropSpeed
Double_t fdSignalPropSpeed
Definition: CbmTofDigitize.h:249
CbmTofDigiBdfPar::GetChanType
Int_t GetChanType(Int_t iSmType, Int_t iRpc) const
Definition: CbmTofDigiBdfPar.cxx:580
CbmDigitizeBase::fCurrentInput
Int_t fCurrentInput
Flag for creation of links to MC.
Definition: CbmDigitizeBase.h:162
CbmTofDigiBdfPar
Parameters class for the CBM ToF digitizer using beam data distributions.
Definition: CbmTofDigiBdfPar.h:30
CbmTofDigitize::fhFiredEvtElCh
TH1 * fhFiredEvtElCh
Definition: CbmTofDigitize.h:328
CbmDigitizeBase::fCurrentEventTime
Double_t fCurrentEventTime
Number of current MC entry.
Definition: CbmDigitizeBase.h:165
CbmTofDigiBdfPar::GetLandauSigma
Double_t GetLandauSigma(Int_t iSmType) const
Definition: CbmTofDigiBdfPar.cxx:655
CbmTofDigitize::fhToTDist
TH1 * fhToTDist
Definition: CbmTofDigitize.h:322
CbmTofDigiBdfPar::GetFeeTimeRes
Double_t GetFeeTimeRes() const
Definition: CbmTofDigiBdfPar.h:49
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmTofDigitize::CreateHistos
Bool_t CreateHistos()
Definition: CbmTofDigitize.cxx:687
CbmTofDetectorId::SetDetectorInfo
virtual Int_t SetDetectorInfo(const CbmTofDetectorInfo detectorInfo)=0
CbmTofCell.h
CbmTofCreateDigiPar::Init
virtual InitStatus Init()
Definition: CbmTofCreateDigiPar.cxx:102
CbmTofDigiPar::GetNode
TGeoNode * GetNode(Int_t iCell)
Definition: CbmTofDigiPar.h:50
first
bool first
Definition: LKFMinuit.cxx:143
CbmTofDigitize::fvlTrckRpcAddr
std::vector< std::vector< ULong64_t > > fvlTrckRpcAddr
Definition: CbmTofDigitize.h:290
CbmTofDigitize::fhMeanFiredPerTrack
TH1 * fhMeanFiredPerTrack
Definition: CbmTofDigitize.h:317
CbmTofDetectorInfo
Definition: CbmTofDetectorId.h:20
CbmTofDigitize::CircleIntersectPosX
Double_t CircleIntersectPosX(Int_t iChanId, Double_t dClustRadius, Double_t dClustCentX, Double_t dClustCentY, Bool_t bUpperSide)
Compute the x position of the intersection of a circle with the upper or.
Definition: CbmTofDigitize.cxx:4685
CbmTofDigitize::fbMcTrkMonitor
Bool_t fbMcTrkMonitor
Definition: CbmTofDigitize.h:346
CbmTofDigitize::LoadBeamtimeValues
Bool_t LoadBeamtimeValues()
Load the beamtime values designed in the parameters: histograms or single values.
Definition: CbmTofDigitize.cxx:401
CbmTofDigitize::DigitizeFlatDisc
Bool_t DigitizeFlatDisc()
Convert TofPoints in input TClonesArray to Tof Digis using an approximation of the.
Definition: CbmTofDigitize.cxx:2084
CbmTofDigitize::fdChannelGain
std::vector< std::vector< std::vector< Double_t > > > fdChannelGain
Definition: CbmTofDigitize.h:260
CbmTofDigitize::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmTofDigitize.cxx:319
CbmTofDigitize::fStorDigi
std::vector< std::vector< std::vector< std::vector< std::pair< CbmTofDigi *, CbmMatch * > > > > > fStorDigi
Definition: CbmTofDigitize.h:280
CbmTofDigitize::fChannelInfo
CbmTofCell * fChannelInfo
Definition: CbmTofDigitize.h:266
CbmTofDigi
Data class for expanded digital TOF information.
Definition: CbmTofDigi.h:38
CbmTofDigitize::fhDigiTimeRes
TH1 * fhDigiTimeRes
Definition: CbmTofDigitize.h:320
CbmTofDigitize::fTimer
TStopwatch fTimer
ROOT timer.
Definition: CbmTofDigitize.h:336
CbmTofDigitize::fStorDigiMatch
std::vector< std::vector< std::vector< std::vector< Int_t > > > > fStorDigiMatch
Definition: CbmTofDigitize.h:282
CbmTofCell::GetSizey
Double_t GetSizey() const
Definition: CbmTofCell.h:41
CbmTofDigiBdfPar::GetLandauMpv
Double_t GetLandauMpv(Int_t iSmType) const
Definition: CbmTofDigiBdfPar.cxx:649
CbmDigitize< CbmTofDigi >::SendData
void SendData(CbmTofDigi *digi, CbmMatch *match=nullptr)
Send a digi and the corresponding match object to the DAQ.
Definition: CbmDigitize.h:219
CbmTofDetectorId::GetCell
virtual Int_t GetCell(const Int_t detectorId)=0
CbmTofDigitize::fhTofPointsPerTrack
TH1 * fhTofPointsPerTrack
Definition: CbmTofDigitize.h:300
CbmTofDigitize::DigitizeDirectClusterSize
Bool_t DigitizeDirectClusterSize()
Convert TofPoints in input TClonesArray to Tof Digis using directly the.
Definition: CbmTofDigitize.cxx:1497
CbmTofGeoHandler::GetCell
Int_t GetCell(Int_t uniqueId)
Definition: CbmTofGeoHandler.cxx:479
CbmTofDigitize::fdNofDigisTot
Double_t fdNofDigisTot
Total number of digis created.
Definition: CbmTofDigitize.h:340
CbmTofDigitize::fTofId
CbmTofDetectorId * fTofId
Definition: CbmTofDigitize.h:264
CbmTofDigitize::fsBeamInputFile
TString fsBeamInputFile
Definition: CbmTofDigitize.h:343
kiNbIntPts
const Int_t kiNbIntPts
Definition: CbmTofDigitize.cxx:59
CbmMCTrack.h
CbmTofDigitize::fhMultiProbElCh
TH1 * fhMultiProbElCh
Definition: CbmTofDigitize.h:329
CbmTofDigitize::fiNofEvents
Int_t fiNofEvents
Total number of events processed.
Definition: CbmTofDigitize.h:337
CbmTofDigitize::fhDigiNbElecCh
TH1 * fhDigiNbElecCh
Definition: CbmTofDigitize.h:314
CbmTofDigitize::DistanceCircleToBase
Double_t DistanceCircleToBase(Double_t dClustRadius, Double_t dBaseXa, Double_t dBaseYa, Double_t dBaseXb, Double_t dBaseYb)
Compute the distance from the cluster center to the base of a disc.
Definition: CbmTofDigitize.cxx:4834
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmTofDetectorId::GetSModule
virtual Int_t GetSModule(const Int_t detectorId)=0
CbmTofDigitize::fsHistoOutFilename
TString fsHistoOutFilename
Definition: CbmTofDigitize.h:298
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTofDigiBdfPar::SetInputFile
void SetInputFile(TString FileName)
Definition: CbmTofDigiBdfPar.h:43
CompTimesExp::data
std::pair< CbmTofDigi *, CbmMatch * > data
Definition: CbmTofDigitize.cxx:64
CbmTofDigiPar
Definition: CbmTofDigiPar.h:18
CbmTofPoint.h
CbmTofGeoHandler::GetSMType
Int_t GetSMType(Int_t uniqueId)
Definition: CbmTofGeoHandler.cxx:459
CbmTofDigiPar::GetCellId
Int_t GetCellId(Int_t i)
Definition: CbmTofDigiPar.h:45
CompTExp
struct CompTimesExp CompTExp
CbmTofDigitize::CompareTimes
Bool_t CompareTimes(CbmTofDigi *p1, CbmTofDigi *p2)
Definition: CbmTofDigitize.cxx:4864
CbmTofDigitize::WriteHistos
Bool_t WriteHistos()
Definition: CbmTofDigitize.cxx:1033
CbmTofDigitize::fhTofPtsInTrkVsGapIndSec
TH2 * fhTofPtsInTrkVsGapIndSec
Definition: CbmTofDigitize.h:303
CbmTofPoint
Geometric intersection of a MC track with a TOFb detector.
Definition: CbmTofPoint.h:40
CbmTofDigitize::fdTimeResElec
Double_t fdTimeResElec
Definition: CbmTofDigitize.h:248
CbmTofDigiBdfPar::GetSignalSpeed
Double_t GetSignalSpeed() const
Definition: CbmTofDigiBdfPar.h:55
CbmTofDigitize::DiscSectionArea
Double_t DiscSectionArea(Double_t dDiscRadius, Double_t dDistBaseToCenter)
Compute area of a disc section from the disc radius and the distance of the.
Definition: CbmTofDigitize.cxx:4661
CbmTofCell::GetY
Double_t GetY() const
Definition: CbmTofCell.h:37
CbmTofDigitize::fvlTrckRpcTime
std::vector< std::vector< Double_t > > fvlTrckRpcTime
Definition: CbmTofDigitize.h:292
CbmTofAddress::GetUniqueAddress
static UInt_t GetUniqueAddress(UInt_t Sm, UInt_t Rpc, UInt_t Channel, UInt_t Side=0, UInt_t SmType=0)
Definition: CbmTofAddress.h:124
CbmTofDigitize::fhDigiMergeTime
TH1 * fhDigiMergeTime
Definition: CbmTofDigitize.h:313
CbmTofDigiBdfPar::GetStartTimeRes
Double_t GetStartTimeRes() const
Definition: CbmTofDigiBdfPar.h:50
CbmTofDigitize::fhNbDigiEvtElCh
TH2 * fhNbDigiEvtElCh
Definition: CbmTofDigitize.h:326
CbmDigitize< CbmTofDigi >::RegisterOutput
void RegisterOutput()
Register the output arrays.
Definition: CbmDigitize.h:175
bFakeBeamCounter
Bool_t bFakeBeamCounter
Definition: CbmTofDigitize.cxx:60
CbmTofDigitize::CircleIntersectPosY
Double_t CircleIntersectPosY(Int_t iChanId, Double_t dClustRadius, Double_t dClustCentX, Double_t dClustCentY, Bool_t bRightSide)
Compute the y position of the intersection of a circle with the left or right edge of a channel....
Definition: CbmTofDigitize.cxx:4758
CbmTofDigiBdfPar::ClearHistos
void ClearHistos()
Definition: CbmTofDigiBdfPar.cxx:79
CbmTofDigitize::fDigiPar
CbmTofDigiPar * fDigiPar
Definition: CbmTofDigitize.h:265
CompTimesExp::operator()
bool operator()(data pair1, data pair2)
Definition: CbmTofDigitize.cxx:65
CbmTofCreateDigiPar
Definition: CbmTofCreateDigiPar.h:23
CbmTofDigitize::CbmTofDigitize
CbmTofDigitize()
Constructor.
Definition: CbmTofDigitize.cxx:85
CbmTofDigiBdfPar::GetClustTotHist
TH1 * GetClustTotHist(Int_t iSmType) const
Definition: CbmTofDigiBdfPar.cxx:643
CbmTofDigitize::fvlTrckChAddr
std::vector< std::vector< ULong64_t > > fvlTrckChAddr
Definition: CbmTofDigitize.h:286
CbmTofDigiBdfPar::GetFeeThreshold
Double_t GetFeeThreshold() const
Definition: CbmTofDigiBdfPar.h:48
CbmTofDigiBdfPar::GetNbSm
Int_t GetNbSm(Int_t iSmType) const
Definition: CbmTofDigiBdfPar.cxx:513
CbmTofDigiBdfPar::LoadBeamtimeHistos
Bool_t LoadBeamtimeHistos()
Definition: CbmTofDigiBdfPar.cxx:314