BmnRoot
Loading...
Searching...
No Matches
BmnGemRaw2Digit.cxx
Go to the documentation of this file.
1#include "BmnGemRaw2Digit.h"
2
3// Boost
4#include <boost/program_options.hpp>
5// BmnRoot
6// #include <BmnEventHeader.h>
7#include <BmnRawDataDecoder.h>
8#include <UniDetectorParameter.h>
9
10namespace po = boost::program_options;
11
12BmnGemRaw2Digit::BmnGemRaw2Digit(Int_t period, Int_t run, TString mapFileName, BmnSetup bmnSetup)
13 : BmnAdcProcessor(period, run, "GEM", ADC_N_CHANNELS, ADC32_N_SAMPLES)
14 , fSmall(nullptr)
15 , fMid(nullptr)
16 , fMapFileName(mapFileName)
17{
18 fSetup = bmnSetup;
19 ReadGlobalMap(fMapFileName);
20 InitArrays();
21
22 fSmall = new BmnGemMap[N_CH_BUF];
23 fMid = new BmnGemMap[N_CH_BUF];
24 Byte_t nMods = 4;
25 for (Byte_t i = 0; i < nMods; i++) {
26 fBigHot.push_back(new BmnGemMap[N_CH_BUF]);
27 fBig.push_back(new BmnGemMap[N_CH_BUF]);
28 }
29
30 if ((fPeriod >= 8) && (fSetup != kSRCSETUP)) {
31 ReadLocalMap("X0_Bottom_Left.txt", fBigHot[2], 2, 2);
32 ReadLocalMap("X0_Bottom_Right.txt", fBigHot[3], 2, 3);
33 ReadLocalMap("X_Bottom_Left.txt", fBig[2], 0, 2);
34 ReadLocalMap("X_Bottom_Right.txt", fBig[3], 0, 3);
35 ReadLocalMap("Y0_Bottom_Left.txt", fBigHot[2], 3, 2);
36 ReadLocalMap("Y0_Bottom_Right.txt", fBigHot[3], 3, 3);
37 ReadLocalMap("Y_Bottom_Left.txt", fBig[2], 1, 2);
38 ReadLocalMap("Y_Bottom_Right.txt", fBig[3], 1, 3);
39
40 ReadLocalMap("X0_Top_Left.txt", fBigHot[0], 2, 0);
41 ReadLocalMap("X0_Top_Right.txt", fBigHot[1], 2, 1);
42 ReadLocalMap("X_Top_Left.txt", fBig[0], 0, 0);
43 ReadLocalMap("X_Top_Right.txt", fBig[1], 0, 1);
44 ReadLocalMap("Y0_Top_Left.txt", fBigHot[0], 3, 0);
45 ReadLocalMap("Y0_Top_Right.txt", fBigHot[1], 3, 1);
46 ReadLocalMap("Y_Top_Left.txt", fBig[0], 1, 0);
47 ReadLocalMap("Y_Top_Right.txt", fBig[1], 1, 1);
48 } else {
49 ReadMap("GEM_X0_middle", fMid, 2, 0);
50 ReadMap("GEM_Y0_middle", fMid, 3, 0);
51 ReadMap("GEM_X1_middle", fMid, 0, 0);
52 ReadMap("GEM_Y1_middle", fMid, 1, 0);
53
54 ReadMap("GEM_X0_Big_Left", fBigHot[1], 2, 1);
55 ReadMap("GEM_Y0_Big_Left", fBigHot[1], 3, 1);
56 ReadMap("GEM_X1_Big_Left", fBig[1], 0, 1);
57 ReadMap("GEM_Y1_Big_Left", fBig[1], 1, 1);
58 ReadMap("GEM_X0_Big_Right", fBigHot[0], 2, 0);
59 ReadMap("GEM_Y0_Big_Right", fBigHot[0], 3, 0);
60 ReadMap("GEM_X1_Big_Right", fBig[0], 0, 0);
61 ReadMap("GEM_Y1_Big_Right", fBig[0], 1, 0);
62 }
63 ReadLocalMap("GEM_X_small.txt", fSmall, 0, 0);
64 ReadLocalMap("GEM_Y_small.txt", fSmall, 1, 0);
65
66 fGemStationSet = BmnGemStripStationSet::Create(period, fSetup);
67 InitNoiseArrays(fGemStationSet);
68}
69
70BmnStatus BmnGemRaw2Digit::ReadGlobalMap(TString FileName)
71{
72 string dummy;
73 Int_t id = 0;
74 UInt_t ser = 0;
75 UInt_t ch_lo = 0;
76 UInt_t ch_hi = 0;
77 UInt_t station = 0;
78 UInt_t mod = 0;
79 UInt_t zone = 0;
80 UInt_t side = 0;
81 UInt_t type = 0;
82
83 vector<string> channel_map;
84 vector<string> spec_thr;
85 Double_t final_thr = 0;
86 Double_t thr_dif = 0;
87 Int_t n_iters = 0;
88 Double_t cmod_cut = 0;
89
90 // Setup options.
91 po::options_description desc("Options");
92 desc.add_options()("CHANNEL_MAP.channels", po::value<vector<string>>(&channel_map)->multitoken(),
93 "ADC channel -> strip map")("SIGNAL_CONFIGURATION.spec_thr",
94 po::value<vector<string>>(&spec_thr)->multitoken(),
95 "ADC channel -> strip map")(
96 "SIGNAL_CONFIGURATION.Threshold", po::value<Double_t>(&final_thr)->default_value(25), "ADC final threshold")(
97 "SIGNAL_CONFIGURATION.ThresholdDif", po::value<Double_t>(&thr_dif)->default_value(10), "ADC final threshold")(
98 "SIGNAL_CONFIGURATION.NIterations", po::value<Int_t>(&n_iters)->default_value(5), "ADC final threshold")(
99 "SIGNAL_CONFIGURATION.CModCut", po::value<Double_t>(&cmod_cut)->default_value(0), "ADC final threshold");
100
101 // Load config file.
102 po::variables_map vm;
103 TString name = TString(getenv("VMCWORKDIR")) + TString("/input/") + fMapFileName;
104 LOG(info) << "Loading GEM Map: Period " << fPeriod << ", Run " << fRun << " " << name;
105 ifstream config_file(name.Data(), ifstream::in);
106 if (!config_file.is_open()) {
107 LOG(error) << "Error opening map-file (" << name << ")!";
108 return kBMNERROR;
109 }
110 po::store(po::parse_config_file(config_file, desc), vm);
111 config_file.close();
112 po::notify(vm);
113
114 // set thresholds
115 SetThreshold(final_thr, thr_dif, n_iters, cmod_cut);
116
117 // set specific thresholds
118 Int_t ch = 0;
119 Double_t thr = 0;
120 for (auto& str : spec_thr) {
121 istringstream ss(str);
122 ss >> std::hex >> ser >> std::dec >> ch >> thr;
123 fSpecThreshMap.insert(make_pair(make_pair(ser, ch), thr));
124 LOGF(debug, "Serial 0x%08X Channel %2d Threshold %3.3f", ser, ch, thr);
125 }
126
127 // set channel map
128 for (auto& str : channel_map) {
129 istringstream ss(str);
130 ss >> std::hex >> ser >> std::dec >> ch_lo >> ch_hi >> id >> side >> type >> station >> mod >> zone;
131 GemMapLine* record = new GemMapLine();
132 record->Ch_hi = ch_hi / GetNSamples();
133 record->Serial = ser;
134 record->Ch_lo = ch_lo / GetNSamples();
135 record->Station = station;
136 record->Zone = zone;
137 record->GEM_id = id;
138 record->Side = side;
139 record->Type = type;
140 record->Module = mod;
141 fMap.push_back(record);
142 GrabNewSerial(ser);
143 }
144 return kBMNSUCCESS;
145}
146
147BmnStatus BmnGemRaw2Digit::ReadLocalMap(TString FileName, BmnGemMap* m, Int_t lay, Int_t mod)
148{
149 TString name = TString(getenv("VMCWORKDIR")) + TString("/input/") + FileName;
150 // printf("%s\n", name.Data());
151 ifstream inFile(name.Data());
152 if (!inFile.is_open()) {
153 LOG(error) << "Error opening map-file (" << name << ")!";
154 return kBMNERROR;
155 }
156 Int_t iStrip = 0;
157 Int_t chan = 0;
158 while (!inFile.eof()) {
159 inFile >> chan;
160 if (!inFile.good())
161 break;
162 m[chan] = BmnGemMap(iStrip++, lay, mod);
163 // printf("\t\t strip %4d lay %4d mod %4d chan %4d\n", iStrip, lay, mod, chan);
164 }
165 return kBMNSUCCESS;
166}
167
168BmnStatus BmnGemRaw2Digit::ReadMap(TString parName, BmnGemMap* m, Int_t lay, Int_t mod)
169{
170 // Int_t size = 0;
171 UniDetectorParameter* par = UniDetectorParameter::GetDetectorParameter("GEM", parName, GetPeriod(), GetRun());
172 vector<UniValue*> iiArr;
173 if (par != NULL)
174 par->GetValue(iiArr);
175 delete par;
176 // printf("%20s mod %d lay %d\n", parName.Data(), mod, lay);
177 for (size_t i = 0; i < iiArr.size(); ++i) {
178 IIValue* pValue = (IIValue*)iiArr[i];
179 m[pValue->value2] = BmnGemMap(pValue->value1 - 1, lay, mod); // Strip begins from 0
180 // printf("\t\t strip %4d lay %4d mod %4d chan %4d %s\n", pValue->value1, lay, mod, pValue->value2,
181 // (pValue->value1 == 0) ? "WFT!!!" : "");
182 }
183
184 if (!iiArr.empty())
185 for (size_t i = 0; i < iiArr.size(); i++)
186 delete iiArr[i];
187 return kBMNSUCCESS;
188}
189
190BmnGemRaw2Digit::~BmnGemRaw2Digit()
191{
192 if (fSmall)
193 delete[] fSmall;
194 if (fMid)
195 delete[] fMid;
196 for (size_t i = 0; i < fBig.size(); i++) {
197 delete[] fBigHot[i];
198 delete[] fBig[i];
199 }
200 if (!fMap.empty())
201 for (size_t i = 0; i < fMap.size(); i++)
202 delete fMap[i];
203 DeleteNoiseArrays(fGemStationSet);
204}
205
206BmnStatus BmnGemRaw2Digit::FillProfiles(TClonesArray* adc)
207{
208 (this->*PrecalcEventModsImp)(adc);
209#ifdef BUILD_DEBUG
210 CalcEventMods();
211#else
212 CalcEventMods_simd();
213#endif
214 ProcessAdc(nullptr, kTRUE);
215 return kBMNSUCCESS;
216}
217
218BmnStatus BmnGemRaw2Digit::FillNoisyChannels()
219{
220 const Int_t kNStripsInBunch = BmnRawDataDecoder::samplesGem;
221 // cout << "GEM: " << kNStripsInBunch << endl;
222 // const Int_t kNBunches = kNStrips / kNStripsInBunch;
223 const double kNThresh = BmnRawDataDecoder::threshGem;
224
225 for (Int_t iCr = 0; iCr < GetNSerials(); ++iCr)
226 for (Int_t iCh = 0; iCh < GetNChannels(); ++iCh)
227 for (auto& it : fMap)
228 if ((Int_t)GetSerials()[iCr] == it->Serial && iCh >= it->Ch_lo && iCh <= it->Ch_hi) {
229 for (Int_t iSmpl = 0; iSmpl < GetNSamples(); ++iSmpl)
230 if (GetNoisyChipChannels()[iCr][iCh][iSmpl] == kTRUE) {
231 Int_t station = -1;
232 Int_t strip = -1;
233 Int_t layer = -1;
234 Int_t mod = -1;
235 MapStrip(it, iCh, iSmpl, station, mod, layer, strip);
236 if (strip < 0)
237 continue;
238 fNoisyChannels[station][mod][layer][strip] = kTRUE;
239 }
240 }
241 for (Int_t iSt = 0; iSt < fGemStationSet->GetNStations(); ++iSt) {
242 if ((iSt == 7) && (fPeriod == 9)) // crutch to avoid marking strips on the beam noisy
243 continue;
244 auto* st = fGemStationSet->GetStation(iSt);
245 for (Int_t iMod = 0; iMod < st->GetNModules(); ++iMod) {
246 auto* mod = st->GetModule(iMod);
247 for (Int_t iLay = 0; iLay < mod->GetNStripLayers(); ++iLay) {
248 TH1F* prof = fSigProf[iSt][iMod][iLay];
249 auto& lay = mod->GetStripLayer(iLay);
250 Int_t kNBunches = lay.GetNStrips() / kNStripsInBunch;
251 for (Int_t iBunch = 0; iBunch < kNBunches; ++iBunch) {
252 Double_t mean = 0.0;
253 for (Int_t iStrip = 0; iStrip < kNStripsInBunch; ++iStrip) {
254 Int_t strip = iStrip + iBunch * kNStripsInBunch;
255 if (fNoisyChannels[iSt][iMod][iLay][strip] == kTRUE)
256 continue;
257 Double_t curr = prof->GetBinContent(strip + 1);
258 // Double_t prev = prof->GetBinContent(strip);
259 mean += curr;
260 }
261 mean /= kNStripsInBunch;
262 for (Int_t iStrip = 0; iStrip < kNStripsInBunch; ++iStrip) {
263 Int_t strip = iStrip + iBunch * kNStripsInBunch;
264 if (fNoisyChannels[iSt][iMod][iLay][strip] == kTRUE)
265 continue;
266 Double_t curr = prof->GetBinContent(strip + 1);
267 // Double_t prev = prof->GetBinContent(strip);
268 // if (kNThresh * meanDiff < curr - prev)
269 if (kNThresh * mean < Abs(curr - mean))
270 fNoisyChannels[iSt][iMod][iLay][strip] = kTRUE;
271 }
272 }
273 }
274 }
275 }
276 for (Int_t iCr = 0; iCr < GetNSerials(); ++iCr)
277 for (Int_t iCh = 0; iCh < GetNChannels(); ++iCh)
278 for (auto& it : fMap)
279 if ((Int_t)GetSerials()[iCr] == it->Serial && iCh >= it->Ch_lo && iCh <= it->Ch_hi)
280 for (Int_t iSmpl = 0; iSmpl < GetNSamples(); ++iSmpl) {
281 Int_t station = -1;
282 Int_t strip = -1;
283 Int_t lay = -1;
284 Int_t mod = -1;
285 MapStrip(it, iCh, iSmpl, station, mod, lay, strip);
286 if (strip < 0)
287 continue;
288 if (fNoisyChannels[station][mod][lay][strip] == kTRUE)
289 GetNoisyChipChannels()[iCr][iCh][iSmpl] = kTRUE;
290 }
291
292 return kBMNSUCCESS;
293}
294
295void BmnGemRaw2Digit::ProcessAdc(TClonesArray* gem, Bool_t doFill)
296{
297 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
298 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
299 for (auto& it : fMap) {
300 // it->Print();
301 if ((Int_t)GetSerials()[iCr] == it->Serial && iCh >= it->Ch_lo && iCh <= it->Ch_hi) {
302 for (Int_t iSmpl = 0; iSmpl < GetNSamples(); ++iSmpl) {
303 if (GetNoisyChipChannels()[iCr][iCh][iSmpl] == kTRUE)
304 continue;
305
306 Int_t station = -1;
307 Int_t strip = -1;
308 Int_t layer = -1;
309 Int_t mod = -1;
310 MapStrip(it, iCh, iSmpl, station, mod, layer, strip);
311 // cout << "DECODER INFO: " << station << " " << mod << " " << layer << " " << strip << endl;
312 if (strip < 0)
313 continue;
314
315 Double_t sig =
316 fAdc[iCr][iCh][iSmpl] - fPedVal[iCr][iCh][iSmpl] + fCMode[iCr][iCh] - fSMode[iCr][iCh];
317 // Double_t Asig = TMath::Abs(sig);
318 Double_t SpecThr(0);
319 auto itThr = fSpecThreshMap.find(make_pair(GetSerials()[iCr], iCh));
320 if (itThr != fSpecThreshMap.end())
321 SpecThr = itThr->second;
322
323 Double_t thr = Max(FinalThr, 3.5 * GetPedestalsRMS()[iCr][iCh][iSmpl]);
324 thr = Max(thr, SpecThr);
325 if ((!GetApplyThreshold()) || (sig > thr)) {
326 if (doFill /* && (Abs(fCMode[iCr][iCh] - fSMode[iCr][iCh]) < cmodcut)*/) {
327 fSigProf[station][mod][layer]->Fill(strip);
328 } else {
329 BmnGemStripDigit* resDig = new ((*gem)[gem->GetEntriesFast()]) BmnGemStripDigit(
330 station, mod, layer, strip, sig, GetPedestalsRMS()[iCr][iCh][iSmpl]);
331 // if ((Abs(fCMode[iCr][iCh] - fSMode[iCr][iCh]) >
332 // cmodcut))
333 // resDig->SetIsGoodDigit(kFALSE);
334 // else
335 resDig->SetIsGoodDigit(kTRUE);
336 }
337 }
338 }
339 break;
340 }
341 }
342 }
343 }
344}
345
346BmnStatus BmnGemRaw2Digit::FillEvent(TClonesArray* adc, TClonesArray* gem)
347{
348 (this->*PrecalcEventModsImp)(adc);
349#ifdef BUILD_DEBUG
350 CalcEventMods();
351#else
352 CalcEventMods_simd();
353#endif
354 ProcessAdc(gem, kFALSE);
355 return kBMNSUCCESS;
356}
357
358inline void BmnGemRaw2Digit::MapStrip(GemMapLine* gemM,
359 UInt_t ch,
360 Int_t iSmpl,
361 Int_t& station,
362 Int_t& mod,
363 Int_t& lay,
364 Int_t& strip)
365{
366 Int_t ch2048 = ch * GetNSamples() + iSmpl;
367 UInt_t realChannel = ch2048;
368 BmnGemMap* fBigMap = nullptr;
369 mod = gemM->Module;
370 UInt_t side = gemM->Side;
371 if (mod > 1)
372 side += 2;
373 if (gemM->GEM_id < 5) {
374 if (gemM->GEM_id < 0) {
375 Int_t chShift = -512;
376 realChannel += chShift;
377 fBigMap = fSmall;
378 } else {
379 fBigMap = fMid;
380 if ((gemM->Ch_hi - gemM->Ch_lo) < 128 / GetNSamples())
381 realChannel = 2048 + ch2048 - gemM->Ch_lo * GetNSamples();
382 }
383 } else {
384 if (gemM->Zone == 0) { // hot zone
385 Int_t chShift = (gemM->Ch_lo == 0) ? 0 : -1024;
386 if (fPeriod < 8)
387 chShift += 1024;
388 // For 6-7 run
389 // Int_t chShift = (gemM->Ch_lo == 0) ? 1024 : 0;
390 realChannel += chShift;
391 fBigMap = fBigHot[side];
392 } else { // big zone
393 fBigMap = fBig[side];
394 if ((gemM->Ch_hi - gemM->Ch_lo) < 128 / GetNSamples())
395 realChannel = 2048 + ch2048 - gemM->Ch_lo * GetNSamples();
396 }
397 }
398 station = gemM->Station;
399 lay = fBigMap[realChannel].lay;
400 strip = fBigMap[realChannel].strip;
401 return;
402}
thr
const float thr
Definition BmnMonProf.cxx:49
ReadMap
int ReadMap(string name)
Definition BmnMonProf.cxx:65
MapStrip
int MapStrip(size_t &iChar, size_t &iCh)
Definition BmnMonProf.cxx:234
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
BmnStatus
BmnStatus
Definition BmnEnums.h:24
kBMNERROR
@ kBMNERROR
Definition BmnEnums.h:26
kBMNSUCCESS
@ kBMNSUCCESS
Definition BmnEnums.h:25
BmnSetup
BmnSetup
Definition BmnEnums.h:89
kSRCSETUP
@ kSRCSETUP
Definition BmnEnums.h:91
BmnAdcProcessor::fSpecThreshMap
map< PlMapKey, Double_t > fSpecThreshMap
Definition BmnAdcProcessor.h:148
BmnAdcProcessor::GetPedestalsRMS
Double_t *** GetPedestalsRMS()
Definition BmnAdcProcessor.h:97
BmnAdcProcessor::GetSerials
vector< UInt_t > & GetSerials()
Definition BmnAdcProcessor.h:109
BmnAdcProcessor::SetThreshold
void SetThreshold(Double_t final_thr, Double_t thr_dif=-1, Double_t n_iters=-1, Double_t cmod_cut=-1)
Definition BmnAdcProcessor.cxx:1207
BmnAdcProcessor::fPedVal
Float_t *** fPedVal
Definition BmnAdcProcessor.h:235
BmnAdcProcessor::fNoisyChannels
Bool_t **** fNoisyChannels
Definition BmnAdcProcessor.h:230
BmnAdcProcessor::GrabNewSerial
void GrabNewSerial(UInt_t serial)
Definition BmnAdcProcessor.cxx:68
BmnAdcProcessor::InitNoiseArrays
void InitNoiseArrays(SST &ss)
Definition BmnAdcProcessor.h:169
BmnAdcProcessor::GetNoisyChipChannels
Bool_t *** GetNoisyChipChannels()
Definition BmnAdcProcessor.h:99
BmnAdcProcessor::DeleteNoiseArrays
void DeleteNoiseArrays(SST &ss)
Definition BmnAdcProcessor.h:202
BmnAdcProcessor::PrecalcEventModsImp
void(BmnAdcProcessor::* PrecalcEventModsImp)(TClonesArray *adc)
Definition BmnAdcProcessor.h:68
BmnGemRaw2Digit::FillEvent
BmnStatus FillEvent(TClonesArray *adc, TClonesArray *gem)
Definition BmnGemRaw2Digit.cxx:346
BmnGemRaw2Digit::BmnGemRaw2Digit
BmnGemRaw2Digit(Int_t period, Int_t run, TString mapFileName, BmnSetup bmnSetup=kBMNSETUP)
Definition BmnGemRaw2Digit.cxx:12
BmnGemRaw2Digit::FillNoisyChannels
BmnStatus FillNoisyChannels()
Definition BmnGemRaw2Digit.cxx:218
BmnGemRaw2Digit::FillProfiles
BmnStatus FillProfiles(TClonesArray *adc)
Definition BmnGemRaw2Digit.cxx:206
BmnGemStripStationSet::Create
static unique_ptr< BmnGemStripStationSet > Create(Int_t period, Int_t stp=0)
Definition BmnGemStripStationSet.cxx:362
BmnRawDataDecoder::threshGem
static double threshGem
Definition BmnRawDataDecoder.h:61
BmnStripDigit::SetIsGoodDigit
void SetIsGoodDigit(Bool_t tmp)
Definition BmnStripDigit.h:40
UniDetectorParameter::GetValue
int GetValue(vector< UniValue * > &parameter_value)
get value of detector parameter presented by an array
Definition UniDetectorParameter.cxx:1107
UniDetectorParameter::GetDetectorParameter
static UniDetectorParameter * GetDetectorParameter(int value_id)
get detector parameter from the database
Definition UniDetectorParameter.cxx:163
ADC_N_CHANNELS
#define ADC_N_CHANNELS
Definition BmnAdcProcessor.h:44
ADC32_N_SAMPLES
#define ADC32_N_SAMPLES
Definition BmnAdcProcessor.h:48
N_CH_BUF
#define N_CH_BUF
Definition BmnGemRaw2Digit.h:25
centrality_clusterization.type
type
Definition centrality_clusterization.py:210
run
-clang-format
setup.name
name
Definition setup.py:7