CbmRoot
Loading...
Searching...
No Matches
CbmPsdSimpleDigitizer.cxx
Go to the documentation of this file.
1/* Copyright (C) 2012-2020 Institute for Nuclear Research, Moscow
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Alla Maevskaya, Selim Seddiki, Sergey Morozov [committer], Volker Friese, Evgeny Kashirin */
4
5// -------------------------------------------------------------------------
6// ----- CbmPsdSimpleDigitizer source file -----
7// ----- Created 15/05/12 by Alla & SELIM & FLORIAN -----
8// ----- Modified 17/03/18 by Sergey Morozov -----
9// -------------------------------------------------------------------------
11
12#include "CbmLink.h"
13#include "CbmMatch.h"
14#include "CbmPsdDigi.h"
15#include "CbmPsdPoint.h"
16
17#include "FairRootManager.h"
18#include <Logger.h>
19
20#include "TClonesArray.h"
21#include "TMath.h"
22#include "TStopwatch.h"
23
24#include <cassert>
25#include <iomanip>
26#include <iostream>
27#include <map>
28
29using std::cout;
30using std::endl;
31using std::fixed;
32using std::right;
33using std::setprecision;
34using std::setw;
35
36using std::map;
37using std::pair;
38
39
40// ----- Default constructor -------------------------------------------
42 : CbmDigitize<CbmPsdDigi>("PsdDigitize")
43 , fNofEvents(0)
44 , fNofPoints(0.)
45 , fNofDigis(0.)
46 , fTimeTot(0.)
47 , fPointArray(NULL)
48{
49}
50// -------------------------------------------------------------------------
51
52
53// ----- Destructor ----------------------------------------------------
55// -------------------------------------------------------------------------
56
57
58// ----- Public method Init --------------------------------------------
60{
61
62 // Get RootManager
63 FairRootManager* ioman = FairRootManager::Instance();
64 assert(ioman),
65
66 // Get input array
67 fPointArray = (TClonesArray*) ioman->GetObject("PsdPoint");
68 assert(fPointArray);
69
70 // Create and register output array
72
73 // --- Read list of inactive channels
74 if (!fInactiveChannelFileName.IsNull()) {
75 LOG(info) << GetName() << ": Reading inactive channels from " << fInactiveChannelFileName;
76 auto result = ReadInactiveChannels();
77 if (!result.second) {
78 LOG(error) << GetName() << ": Error in reading from file! Task will be inactive.";
79 return kFATAL;
80 }
81 LOG(info) << GetName() << ": " << std::get<0>(result) << " lines read from file, " << fInactiveChannels.size()
82 << " channels set inactive";
83 }
84
85 // Statistics
86 fNofEvents = 0;
87 fNofPoints = 0;
88 fNofDigis = 0.;
89 fTimeTot = 0.;
90
91 LOG(info) << fName << ": Initialisation successful " << kSUCCESS;
92 return kSUCCESS;
93}
94// -------------------------------------------------------------------------
95
96
97// ----- Public method Exec --------------------------------------------
99{
100 // Event info (for event time)
101 GetEventInfo();
102
103 TStopwatch timer;
104 timer.Start();
105
106 // Loop over PsdPoints
107 Int_t nPoints = fPointArray->GetEntriesFast();
108 LOG(debug) << fName << ": processing event " << fCurrentEvent << " at t = " << fCurrentEventTime << " ns";
109
110 // Declare some variables
111 CbmPsdPoint* point = nullptr;
112 Int_t modID = -1; // module ID
113 Int_t scinID = -1; // #sciillator
114 Int_t sec;
115 std::map<UInt_t, std::pair<CbmPsdDigi, CbmMatch*>> fired_digis_map; //store Digis and Matches for each module
116
117 for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
118 point = (CbmPsdPoint*) fPointArray->At(iPoint);
119 if (!point) continue;
120
121 CbmLink link(1., iPoint, fCurrentMCEntry, fCurrentInput);
122
123 modID = point->GetModuleID(); //marina 1-44 (45)
124 scinID = point->GetDetectorID(); //1-60
125 Double_t eLoss = point->GetEnergyLoss();
126 Double_t pTime = point->GetTime();
127 sec = (Int_t)((scinID - 1) / 6) + 1; //marina 1-10
128 UInt_t uAddress = CbmPsdAddress::GetAddress(modID, sec);
129
130 auto it = fired_digis_map.find(uAddress);
131 if (it != fired_digis_map.end()) {
132 //this key exists
133 it->second.first.SetEdep(it->second.first.GetEdep() + eLoss);
134 if (pTime < it->second.first.GetTime()) it->second.first.SetTime(pTime);
135 it->second.second->AddLink(link);
136 }
137 else {
138 //this key is new
139 CbmPsdDigi digi = CbmPsdDigi(uAddress, pTime, eLoss);
140 CbmMatch* match = new CbmMatch();
141 match->AddLink(link);
142 fired_digis_map.insert(std::make_pair(uAddress, std::make_pair(digi, match)));
143 }
144 } // Loop over MCPoints
145
146 Int_t nDigis = 0;
147 for (auto entry : fired_digis_map) {
148 Double_t eDep = entry.second.first.GetEdep();
149 Double_t eLossMIP = eDep / 0.005; // 5MeV per MIP
150 Double_t pixPerMIP = 15.; // 15 pix per MIP
151 Double_t eLossMIPSmeared = gRandom->Gaus(eLossMIP * pixPerMIP, sqrt(eLossMIP * pixPerMIP)) / pixPerMIP;
152 Double_t eLossSmeared = eLossMIPSmeared * 0.005;
153 Double_t eNoise = gRandom->Gaus(0, 15) / 50. * 0.005;
154 eLossSmeared += eNoise;
155 CbmPsdDigi* digi =
156 new CbmPsdDigi(entry.second.first.GetAddress(), entry.second.first.GetTime() + fCurrentEventTime, eLossSmeared);
157 SendData(digi->GetTime(), digi, entry.second.second);
158 nDigis++;
159 LOG(debug1) << fName << ": Digi " << nDigis << " Time " << entry.second.first.GetTime() + fCurrentEventTime
160 << " Section " << entry.second.first.GetSectionID() << " Module " << entry.second.first.GetModuleID()
161 << " energy " << eLossSmeared;
162 }
163
164 // --- Event log
165 timer.Stop();
166 LOG(info) << "+ " << setw(15) << GetName() << ": Event " << setw(6) << right << fCurrentEvent << " at " << fixed
167 << setprecision(3) << fCurrentEventTime << " ns, points: " << nPoints << ", digis: " << nDigis
168 << ". Exec time " << setprecision(6) << timer.RealTime() << " s.";
169
170 // --- Run statistics
171 fNofEvents++;
172 fNofPoints += nPoints;
173 fNofDigis += nDigis;
174 fTimeTot += timer.RealTime();
175}
176// -------------------------------------------------------------------------
177
178
179// ----- End-of-run ----------------------------------------------------
181{
182 std::cout << std::endl;
183 LOG(info) << "=====================================";
184 LOG(info) << GetName() << ": Run summary";
185 LOG(info) << "Events processed : " << fNofEvents;
186 LOG(info) << "PsdPoint / event : " << setprecision(1) << fNofPoints / Double_t(fNofEvents);
187 LOG(info) << "PsdDigi / event : " << fNofDigis / Double_t(fNofEvents);
188 LOG(info) << "Digis per point : " << setprecision(6) << fNofDigis / fNofPoints;
189 LOG(info) << "Real time per event : " << fTimeTot / Double_t(fNofEvents) << " s";
190 LOG(info) << "=====================================";
191}
192// -------------------------------------------------------------------------
193
194
ClassImp(CbmConverterManager)
friend fvec sqrt(const fvec &a)
virtual std::pair< size_t, bool > ReadInactiveChannels()
Set of inactive channels, indicated by CbmAddress.
Int_t fCurrentEvent
Number of current input.
TString fInactiveChannelFileName
Time of current MC event [ns].
std::set< uint32_t > fInactiveChannels
Name of file with inactive channels.
void GetEventInfo()
Get event information.
Int_t fCurrentInput
Start time of run [ns].
Double_t fCurrentEventTime
Number of current MC entry.
Int_t fCurrentMCEntry
Number of current MC event.
Base class template for CBM digitisation tasks.
Definition CbmDigitize.h:44
void SendData(Double_t time, CbmPsdDigi *digi, CbmMatch *match=nullptr)
void AddLink(const CbmLink &newLink)
Definition CbmMatch.cxx:47
static uint32_t GetAddress(int32_t moduleId, int32_t sectionId)
Return address from system ID, module, Section.
Data class for PSD digital information.
Definition CbmPsdDigi.h:36
double GetTime() const
Time.
Definition CbmPsdDigi.h:105
int32_t GetModuleID() const
Definition CbmPsdPoint.h:61
virtual void Exec(Option_t *opt)
virtual void Finish()
End-of-run action.