SND@LHC Software
Loading...
Searching...
No Matches
ConvRawData.cxx
Go to the documentation of this file.
1#include <iostream>
2#include <algorithm>
3#include <fstream>
4#include <sstream>
5#include <map>
6#include <tuple> // to return multiple values
7#include <chrono> // for time estimations
8#include <ctime>
9#include <string>
10#include <stdint.h>
11#include <stdlib.h> // exit
12#include <vector>
13#include <algorithm> // std::sort
14#include <sys/stat.h>
15#include "TKey.h"
16#include "FairEventHeader.h" // for FairEventHeader
17#include "SNDLHCEventHeader.h" // for EventHeader
18#include "FairRunAna.h" // for FairRunAna
19#include "FairRootManager.h" // for FairRootManager
20#include "FairParamList.h"
21#include "FairLogger.h"
22#include "ConvRawData.h"
23#include "TClonesArray.h"
24#include "TObjString.h"
25#include "TROOT.h"
26#include "TList.h"
27#include "TFile.h"
28#include "TTree.h"
29#include "TBranch.h"
30#include "TLeaf.h"
31#include "TH1F.h"
32#include "TStopwatch.h"
33#include "TSystem.h"
34#include "SndlhcHit.h"
35#include "Scifi.h" // for SciFi detector
36#include "sndScifiHit.h" // for SciFi Hit
37#include "MuFilterHit.h" // for Muon Filter Hit
38#include "TString.h"
39#include "nlohmann/json.hpp" // library to operate with json files
40#include "boardMappingParser.h" // for board mapping
41#include "XrdCl/XrdClFile.hh"
42
43using namespace std;
44using namespace std::chrono;
45using namespace XrdCl;
46
47
48ConvRawData::ConvRawData()
49 : FairTask("ConvRawData")
50 , fEventTree(nullptr)
51 , boards{}
52 , fSNDLHCEventHeader(nullptr)
53 , fEventHeader(nullptr)
54 , fDigiSciFi(nullptr)
55 , fDigiMuFilter(nullptr)
56{}
57
58ConvRawData::~ConvRawData() {}
59
60InitStatus ConvRawData::Init()
61{
62 FairRootManager* ioman = FairRootManager::Instance();
63 if (!ioman) {
64 LOG (error) << "-E- ConvRawData::Init: "
65 << "RootManager not instantiated!";
66 return kFATAL;
67 }
68
69 // Reading input - initiating parameters
70 TObjString* runN_obj = dynamic_cast<TObjString*>(ioman->GetObject("runN"));
71 TObjString* pathCalib_obj = dynamic_cast<TObjString*>(ioman->GetObject("pathCalib"));
72 TObjString* pathJSON_obj = dynamic_cast<TObjString*>(ioman->GetObject("pathJSON"));
73 TObjString* nEvents_obj = dynamic_cast<TObjString*>(ioman->GetObject("nEvents"));
74 TObjString* nStart_obj = dynamic_cast<TObjString*>(ioman->GetObject("nStart"));
75 TObjString* debug_obj = dynamic_cast<TObjString*>(ioman->GetObject("debug"));
76 TObjString* stop_obj = dynamic_cast<TObjString*>(ioman->GetObject("stop"));
77 TObjString* heartBeat_obj = dynamic_cast<TObjString*>(ioman->GetObject("heartBeat"));
78 TObjString* makeCalibration_obj = dynamic_cast<TObjString*>(ioman->GetObject("makeCalibration"));
79 TObjString* chi2Max_obj = dynamic_cast<TObjString*>(ioman->GetObject("chi2Max"));
80 TObjString* saturationLimit_obj = dynamic_cast<TObjString*>(ioman->GetObject("saturationLimit"));
81 TObjString* newFormat_obj = dynamic_cast<TObjString*>(ioman->GetObject("newFormat"));
82 TObjString* local_obj = dynamic_cast<TObjString*>(ioman->GetObject("local"));
83 // TMap containing the filling scheme per given run
84 FSmap = static_cast<TMap*>(ioman->GetObject("FSmap"));
85 // Input raw data file is read from the FairRootManager
86 // This allows to have it in custom format, e.g. have arbitary names of TTrees
87 TFile* f0 = dynamic_cast<TFile*>(ioman->GetObject("rawData"));
88 // Set run parameters
89 std::istringstream(runN_obj->GetString().Data()) >> frunNumber;
90 fpathCalib = pathCalib_obj->GetString().Data();
91 fpathJSON = pathJSON_obj->GetString().Data();
92 std::istringstream(nEvents_obj->GetString().Data()) >> fnEvents;
93 std::istringstream(nStart_obj->GetString().Data()) >> fnStart;
94 std::istringstream(debug_obj->GetString().Data()) >> debug;
95 std::istringstream(stop_obj->GetString().Data()) >> stop;
96 std::istringstream(heartBeat_obj->GetString().Data()) >> fheartBeat;
97 std::istringstream(makeCalibration_obj->GetString().Data()) >> makeCalibration;
98 std::istringstream(chi2Max_obj->GetString().Data()) >> chi2Max;
99 std::istringstream(saturationLimit_obj->GetString().Data()) >> saturationLimit;
100 std::istringstream(newFormat_obj->GetString().Data()) >> newFormat;
101 std::istringstream(local_obj->GetString().Data()) >> local;
102
103 if (!newFormat)
104 {
105 fEventTree = (TTree*)f0->Get("event");
106 // Get board_x data
107 TIter next(f0->GetListOfKeys());
108 TKey *b = new TKey();
109 string name;
110 while ((b = (TKey*)next()))
111 {
112 name = b->GetName();
113 // string.find func: If no matches were found, the function returns string::npos.
114 if ( name.find("board") == string::npos ) continue;
115 boards[name] = (TTree*)f0->Get(name.c_str());
116 }
117 // use FairRoot eventHeader class
118 fEventHeader = new FairEventHeader();
119 ioman->Register("EventHeader", "sndEventHeader", fEventHeader, kTRUE);
120 }
121 else
122 {
123 fEventTree = (TTree*)f0->Get("data");
124 // use sndlhc eventHeader class
125 fSNDLHCEventHeader = new SNDLHCEventHeader();
126 ioman->Register("EventHeader.", "sndEventHeader", fSNDLHCEventHeader, kTRUE);
127 }
128
129 fDigiSciFi = new TClonesArray("sndScifiHit");
130 ioman->Register("Digi_ScifiHits", "DigiScifiHit_det", fDigiSciFi, kTRUE);
131 fDigiMuFilter = new TClonesArray("MuFilterHit");
132 ioman->Register("Digi_MuFilterHits", "DigiMuFilterHit_det", fDigiMuFilter, kTRUE);
133 ScifiDet = dynamic_cast<Scifi*> (gROOT->GetListOfGlobals()->FindObject("Scifi") );
134 TStopwatch timerCSV;
135 timerCSV.Start();
136 read_csv(fpathCalib);
137 timerCSV.Stop();
138 LOG (info) << "Time to read CSV " << timerCSV.RealTime();
139
140 //calibrationReport();
141 TStopwatch timerBMap;
142 timerBMap.Start();
143 DetMapping(fpathJSON);
144 if (newFormat) StartTimeofRun(fpathJSON);
145 timerBMap.Stop();
146 LOG (info) << "Time to set the board mapping " << timerBMap.RealTime();
147
148 eventNumber = fnStart;
149
150 return kSUCCESS;
151}
152
153void ConvRawData::Exec(Option_t* /*opt*/)
154{
155 fDigiSciFi->Delete();
156 fDigiMuFilter->Delete();
157
158 if (!newFormat) Process0();
159 else Process1();
160
161 // Manually change event number as input file is not set as source for this task
162 eventNumber++;
163
164}
165
166void ConvRawData::UpdateInput(int NewStart)
167{
168 fEventTree->Refresh();
169 if (!newFormat)
170 for (auto it : boards) boards[it.first]->Refresh();
171 eventNumber = NewStart;
172}
173
174void ConvRawData::Process0()
175{
176 int indexSciFi{}, indexMuFilter{};
177 bool scifi = false, mask = false;
178 int board_id{};
179 TTree* bt;
180 int tofpet_id{}, tofpet_channel{}, tac{}, mat{};
181 string station;
182 double TDC{}, QDC{}, Chi2ndof{}, satur{};
183 int A{};
184 double B{};
185 high_resolution_clock::time_point tE{}, t0{}, t1{}, t4{}, t5{}, t6{}, tt{};
186 int system{}, key{}, sipmChannel{};
187 string tmp;
188 int nSiPMs{}, nSides{}, direction{}, detID{}, sipm_number{}, chan{}, orientation{}, sipmLocal{};
189 int sipmID{};
190 double test{};
191 //TStopwatch timer;
192
193 tE = high_resolution_clock::now();
194 //timer.Start();
195 fEventTree->GetEvent(eventNumber);
196 if ( eventNumber%fheartBeat == 0 )
197 {
198 tt = high_resolution_clock::now();
199 time_t ttp = high_resolution_clock::to_time_t(tt);
200 LOG (info) << "run " << frunNumber << " event " << eventNumber
201 << " local time " << ctime(&ttp);
202 }
203
204 fEventHeader->SetRunId(frunNumber);
205 fEventHeader->SetEventTime(fEventTree->GetLeaf("evtTimestamp")->GetValue());
206 LOG (info) << "event: " << eventNumber << " timestamp: "
207 << fEventTree->GetLeaf("evtTimestamp")->GetValue();
208 // Delete pointer map elements
209 for (auto it : digiSciFiStore)
210 {
211 delete it.second;
212 }
213 digiSciFiStore.clear();
214 for (auto it : digiMuFilterStore)
215 {
216 delete it.second;
217 }
218 digiMuFilterStore.clear();
219
220 // Loop over boards
221 for ( auto board : boards )// loop over TTrees
222 {
223 board_id = stoi(board.first.substr(board.first.find("_")+1));
224 scifi = true;
225 if (boardMaps["Scifi"].count(board.first)!=0)
226 {
227 for (auto it : boardMaps["Scifi"][board.first])
228 {
229 station = it.first;
230 mat = it.second;
231 }
232 }
233 else if (boardMaps["MuFilter"].count(board.first)!=0) scifi = false;
234 else
235 {
236 LOG (error) << board.first << " not known. Serious error, stop!";
237 break;
238 }
239 bt = boards[board.first];
240 bt->GetEvent(eventNumber);
241 // Loop over hits in event
242 for ( int n = 0; n < bt->GetLeaf("nHits")->GetValue(); n++ )
243 {
244 mask = false;
245 LOG (info) << "In scifi? " << scifi
246 << " " << board_id << " " << bt->GetLeaf("tofpetId")->GetValue(n)
247 << " " << bt->GetLeaf("tofpetChannel")->GetValue(n)
248 << " " << bt->GetLeaf("tac")->GetValue(n)
249 << " " << bt->GetLeaf("tCoarse")->GetValue(n)
250 << " " << bt->GetLeaf("tFine")->GetValue(n)
251 << " " << bt->GetLeaf("vCoarse")->GetValue(n)
252 << " " << bt->GetLeaf("vFine")->GetValue(n);
253 t0 = high_resolution_clock::now();
254 tofpet_id = bt->GetLeaf("tofpetId")->GetValue(n);
255 tofpet_channel = bt->GetLeaf("tofpetChannel")->GetValue(n);
256 tac = bt->GetLeaf("tac")->GetValue(n);
257 /* Since run 39 calibration is done online and
258 calib data stored in root raw data file */
259 if (makeCalibration)
260 tie(TDC,QDC,Chi2ndof,satur) = comb_calibration(board_id, tofpet_id, tofpet_channel, tac,
261 bt->GetLeaf("vCoarse")->GetValue(n),
262 bt->GetLeaf("vFine")->GetValue(n),
263 bt->GetLeaf("tCoarse")->GetValue(n),
264 bt->GetLeaf("tFine")->GetValue(n),
265 1.0, 0);
266 else
267 {
268 TDC = bt->GetLeaf("timestamp")->GetValue(n);
269 QDC = bt->GetLeaf("value")->GetValue(n);
270 Chi2ndof = max(bt->GetLeaf("timestampCalChi2")->GetValue(n)/bt->GetLeaf("timestampCalDof")->GetValue(n),
271 bt->GetLeaf("valueCalChi2")->GetValue(n)/bt->GetLeaf("valueCalDof")->GetValue(n));
272 satur = bt->GetLeaf("valueSaturation")->GetValue(n);
273 }
274
275 // Print a warning if TDC or QDC is nan.
276 if ( TDC != TDC || QDC!=QDC) {
277 LOG (error) << "NAN tdc/qdc detected! Check maps!"
278 << " " << board_id << " " << bt->GetLeaf("tofpetId")->GetValue(n)
279 << " " << bt->GetLeaf("tofpetChannel")->GetValue(n)
280 << " " << bt->GetLeaf("tac")->GetValue(n)
281 << " " << bt->GetLeaf("tCoarse")->GetValue(n)
282 << " " << bt->GetLeaf("tFine")->GetValue(n)
283 << " " << bt->GetLeaf("vCoarse")->GetValue(n)
284 << " " << bt->GetLeaf("vFine")->GetValue(n);
285 }
286
287 t1 = high_resolution_clock::now();
288 if ( Chi2ndof > chi2Max )
289 {
290 if (QDC>1E20) QDC = 997.; // checking for inf
291 if (QDC != QDC) QDC = 998.; // checking for nan
292 if (QDC>0) QDC = -QDC;
293 mask = true;
294 }
295 else if (satur > saturationLimit || QDC>1E20 || QDC != QDC)
296 {
297 if (QDC>1E20) QDC = 987.; // checking for inf
298 LOG (info) << "inf " << board_id << " " << bt->GetLeaf("tofpetId")->GetValue(n)
299 << " " << bt->GetLeaf("tofpetChannel")->GetValue(n)
300 << " " << bt->GetLeaf("tac")->GetValue(n)
301 << " " << bt->GetLeaf("vCoarse")->GetValue(n)
302 << " " << bt->GetLeaf("vFine")->GetValue(n)
303 << " " << TDC-bt->GetLeaf("tCoarse")->GetValue(n)
304 << " " << eventNumber << " " << Chi2ndof;
305 if (QDC != QDC) QDC = 988.; // checking for nan
306 LOG (info) << "nan " << board_id << " " << bt->GetLeaf("tofpetId")->GetValue(n)
307 << " " << bt->GetLeaf("tofpetChannel")->GetValue(n)
308 << " " << bt->GetLeaf("tac")->GetValue(n)
309 << " " << bt->GetLeaf("vCoarse")->GetValue(n)
310 << " " << bt->GetLeaf("vFine")->GetValue(n)
311 << " " << TDC-bt->GetLeaf("tCoarse")->GetValue(n)
312 << " " << TDC << " " << bt->GetLeaf("tCoarse")->GetValue(n)
313 << " " << eventNumber << " " << Chi2ndof;
314 A = int(min(QDC,double(1000.)));
315 B = min(satur,double(999.))/1000.;
316 QDC = A+B;
317 mask = true;
318 }
319 else if ( Chi2ndof > chi2Max )
320 {
321 if (QDC>0) QDC = -QDC;
322 mask = true;
323 }
324 LOG (info) << "calibrated: tdc = " << TDC << ", qdc = " << QDC;//TDC clock cycle = 160 MHz or 6.25ns
325 t4 = high_resolution_clock::now();
326 // Set the unit of the execution time measurement to ns
327 counters["qdc"]+= duration_cast<nanoseconds>(t1 - t0).count();
328 counters["make"]+= duration_cast<nanoseconds>(t4-t0).count();
329
330 // MuFilter encoding
331 if (!scifi)
332 {
333 system = MufiSystem[board_id][tofpet_id];
334 key = (tofpet_id%2)*1000 + tofpet_channel;
335 tmp = boardMapsMu["MuFilter"][board.first][slots[tofpet_id]];
336 if (debug || !(tmp.find("not") == string::npos))
337 {
338 LOG (info) << system << " " << key << " " << board.first << " " << tofpet_id
339 << " " << tofpet_id%2 << " " << tofpet_channel;
340 }
341 sipmChannel = 99;
342 if ( TofpetMap[system].count(key) == 0)
343 {
344 cout << "key " << key << " does not exist. " << endl
345 << board.first << " Tofpet id " << tofpet_id
346 << " System " << system << " has Tofpet map elements: {";
347 for ( auto it : TofpetMap[system])
348 {
349 cout << it.first << " : " << it.second << ", ";
350 }
351 cout << "}\n";
352 }
353 else sipmChannel = TofpetMap[system][key]-1;
354
355 nSiPMs = abs(offMap[tmp][1]);
356 nSides = abs(offMap[tmp][2]);
357 direction = int(offMap[tmp][1]/nSiPMs);
358 detID = offMap[tmp][0] + direction*int(sipmChannel/nSiPMs);
359 sipm_number = sipmChannel%(nSiPMs);
360 if ( tmp.find("Right") != string::npos ) sipm_number+= nSiPMs;
361 if (digiMuFilterStore.count(detID)==0)
362 {
363 digiMuFilterStore[detID] = new MuFilterHit(detID,nSiPMs,nSides);
364 }
365 test = digiMuFilterStore[detID]->GetSignal(sipm_number);
366 digiMuFilterStore[detID]->SetDigi(QDC,TDC,sipm_number);
367 digiMuFilterStore[detID]->SetDaqID(sipm_number,n, board_id, tofpet_id, tofpet_channel);
368 if (mask) digiMuFilterStore[detID]->SetMasked(sipm_number);
369
370 LOG (info) << "create mu hit: " << detID << " " << tmp << " " << system
371 << " " << tofpet_id << " " << offMap[tmp][0] << " " << offMap[tmp][1]
372 << " " << offMap[tmp][2] << " " << sipmChannel << " " << nSiPMs
373 << " " << nSides << " " << test << endl
374 << detID << " " << sipm_number << " " << QDC << " " << TDC;
375
376 if (test>0 || detID%1000>200 || sipm_number>15)
377 {
378 cout << "what goes wrong? " << detID << " SiPM " << sipm_number << " system " << system
379 << " key " << key << " board " << board.first << " tofperID " << tofpet_id
380 << " tofperChannel " << tofpet_channel << " test " << test << endl;
381 }
382 t5 = high_resolution_clock::now();
383 counters["createMufi"]+= duration_cast<nanoseconds>(t5 - t4).count();
384 } // end MuFilter encoding
385
386 else // now Scifi encoding
387 {
388 chan = channel_func(tofpet_id, tofpet_channel, mat);
389 orientation = 1;
390 if (station[2]=='Y') orientation = 0;
391 sipmLocal = (chan - mat*512);
392 sipmID = 1000000*int(station[1]-'0') + 100000*orientation + 10000*mat
393 + 1000*(int(sipmLocal/128)) + chan%128;
394 if (digiSciFiStore.count(sipmID)==0)
395 {
396 digiSciFiStore[sipmID] = new sndScifiHit(sipmID);
397 }
398 digiSciFiStore[sipmID]->SetDigi(QDC,TDC);
399 digiSciFiStore[sipmID]->SetDaqID(0,n, board_id, tofpet_id, tofpet_channel);
400 if (mask) digiSciFiStore[sipmID]->setInvalid();
401 LOG (info) << "create scifi hit: tdc = " << board.first << " " << sipmID
402 << " " << QDC << " " << TDC <<endl
403 << "tofpet:" << " " << tofpet_id << " " << tofpet_channel << " " << mat
404 << " " << chan << endl
405 << station[1] << " " << station[2] << " " << mat << " " << chan
406 << " " << int(chan/128)%4 << " " << chan%128;
407 t5 = high_resolution_clock::now();
408 counters["createScifi"]+= duration_cast<nanoseconds>(t5 - t4).count();
409 } // end Scifi encoding
410 } // end loop over hits in event
411 } // end loop over boards (TTrees in data file)
412
413 counters["N"]+= 1;
414 t6 = high_resolution_clock::now();
415 for (auto it_sipmID : digiSciFiStore)
416 {
417 (*fDigiSciFi)[indexSciFi]=digiSciFiStore[it_sipmID.first];
418 indexSciFi+= 1;
419 }
420 for (auto it_detID : digiMuFilterStore)
421 {
422 (*fDigiMuFilter)[indexMuFilter]=digiMuFilterStore[it_detID.first];
423 indexMuFilter+= 1;
424 }
425 counters["storage"]+= duration_cast<nanoseconds>(high_resolution_clock::now() - t6).count();
426 counters["event"]+= duration_cast<nanoseconds>(high_resolution_clock::now() - tE).count();
427 //timer.Stop();
428 //cout<<timer.RealTime()<<endl;
429
430 LOG (info) << fnStart+1 << " events processed out of "
431 << fEventTree->GetEntries() << " number of events in file.";
432 /*
433 cout << "Monitor computing time:" << endl;
434 for (auto it: counters)
435 {
436 if( it.first=="N" )
437 {
438 cout << "Processed " << it.first<< " = " << it.second << " events." << endl;
439 }
440 else
441 {
442 // Print execution time in seconds
443 cout << "stage " << it.first<< " took " << it.second*1e-9 << " [s]" << endl;
444 }
445 }
446 */
447}
448
449void ConvRawData::Process1()
450{
451 int indexSciFi{}, indexMuFilter{};
452 bool scifi = false, mask = false;
453 int board_id{};
454 string board_name;
455 int tofpet_id{}, tofpet_channel{}, tac{}, mat{};
456 string station;
457 double TDC{}, QDC{}, Chi2ndof{}, satur{};
458 int A{};
459 double B{};
460 high_resolution_clock::time_point tE{}, t0{}, t1{}, t4{}, t5{}, t6{}, tt{};
461 int system{}, key{}, sipmChannel{};
462 string tmp;
463 int nSiPMs{}, nSides{}, direction{}, detID{}, sipm_number{}, chan{}, orientation{}, sipmLocal{};
464 int sipmID{};
465 double test{};
466 int64_t eventTime{};
467 //TStopwatch timer;
468
469 tE = high_resolution_clock::now();
470 //timer.Start();
471 fEventTree->GetEvent(eventNumber);
472 if ( eventNumber%fheartBeat == 0 )
473 {
474 tt = high_resolution_clock::now();
475 time_t ttp = high_resolution_clock::to_time_t(tt);
476 LOG (info) << "run " << frunNumber << " event " << eventNumber
477 << " local time " << ctime(&ttp);
478 }
479
480 fSNDLHCEventHeader->SetFlags(fEventTree->GetLeaf("evtFlags")->GetValue());
481 fSNDLHCEventHeader->SetRunId(frunNumber);
482 eventTime = fEventTree->GetLeaf("evtTimestamp")->GetValue();
483 fSNDLHCEventHeader->SetEventTime(eventTime);
484 fSNDLHCEventHeader->SetUTCtimestamp(eventTime*6.23768*1e-9 + runStartUTC);
485 fSNDLHCEventHeader->SetEventNumber(fEventTree->GetLeaf("evtNumber")->GetValue());
486 // Fill filling scheme data into the event header
487 if (FSmap->GetEntries()>1){
488 // ion runs
489 if (fSNDLHCEventHeader->GetAccMode()== 12){
490 fSNDLHCEventHeader->SetBunchType(stoi(((TObjString*)FSmap->GetValue(Form("%d",
491 int((eventTime%(4*3564))/8+0.5))))->GetString().Data()));
492 }
493 // proton runs
494 else{
495 fSNDLHCEventHeader->SetBunchType(stoi(((TObjString*)FSmap->GetValue(Form("%d",
496 int((eventTime%(4*3564))/4))))->GetString().Data()));
497 }
498 }
499 else
500 fSNDLHCEventHeader->SetBunchType(stoi(((TObjString*)FSmap->GetValue("0"))->GetString().Data()));
501
502 LOG (info) << "evtNumber per run "
503 << fEventTree->GetLeaf("evtNumber")->GetValue()
504 << " evtNumber per partition: " << eventNumber
505 << " timestamp: " << eventTime;
506 // Delete pointer map elements
507 for (auto it : digiSciFiStore)
508 {
509 delete it.second;
510 }
511 digiSciFiStore.clear();
512 for (auto it : digiMuFilterStore)
513 {
514 delete it.second;
515 }
516 digiMuFilterStore.clear();
517
518 // Loop over hits per event!
519 for ( int n =0; n < fEventTree->GetLeaf("nHits")->GetValue(); n++ )
520 {
521 board_id = fEventTree->GetLeaf("boardId")->GetValue(n);
522 board_name = "board_"+to_string(board_id);
523 scifi = true;
524 if (boardMaps["Scifi"].count(board_name)!=0)
525 {
526 for (auto it : boardMaps["Scifi"][board_name])
527 {
528 station = it.first;
529 mat = it.second;
530 }
531 }
532 else if (boardMaps["MuFilter"].count(board_name)!=0) scifi = false;
533 else
534 {
535 LOG (error) << board_name << " not known. Serious error, stop!";
536 break;
537 }
538 mask = false;
539 LOG (info) << "In scifi? " << scifi
540 << " " << board_id << " " << fEventTree->GetLeaf("tofpetId")->GetValue(n)
541 << " " << fEventTree->GetLeaf("tofpetChannel")->GetValue(n)
542 << " " << fEventTree->GetLeaf("tac")->GetValue(n)
543 << " " << fEventTree->GetLeaf("tCoarse")->GetValue(n)
544 << " " << fEventTree->GetLeaf("tFine")->GetValue(n)
545 << " " << fEventTree->GetLeaf("vCoarse")->GetValue(n)
546 << " " << fEventTree->GetLeaf("vFine")->GetValue(n);
547 t0 = high_resolution_clock::now();
548 tofpet_id = fEventTree->GetLeaf("tofpetId")->GetValue(n);
549 tofpet_channel = fEventTree->GetLeaf("tofpetChannel")->GetValue(n);
550 tac = fEventTree->GetLeaf("tac")->GetValue(n);
551 /* Since run 39 calibration is done online and
552 calib data stored in root raw data file */
553 if (makeCalibration)
554 tie(TDC,QDC,Chi2ndof,satur) = comb_calibration(board_id, tofpet_id, tofpet_channel, tac,
555 fEventTree->GetLeaf("vCoarse")->GetValue(n),
556 fEventTree->GetLeaf("vFine")->GetValue(n),
557 fEventTree->GetLeaf("tCoarse")->GetValue(n),
558 fEventTree->GetLeaf("tFine")->GetValue(n),
559 1.0, 0);
560 else
561 {
562 TDC = fEventTree->GetLeaf("timestamp")->GetValue(n);
563 QDC = fEventTree->GetLeaf("value")->GetValue(n);
564 Chi2ndof = max(fEventTree->GetLeaf("timestampCalChi2")->GetValue(n)/fEventTree->GetLeaf("timestampCalDof")->GetValue(n),
565 fEventTree->GetLeaf("valueCalChi2")->GetValue(n)/fEventTree->GetLeaf("valueCalDof")->GetValue(n));
566 satur = fEventTree->GetLeaf("valueSaturation")->GetValue(n);
567 }
568
569 // Print a warning if TDC or QDC is nan.
570 if ( TDC != TDC || QDC!=QDC) {
571 LOG (error) << "NAN tdc/qdc detected! Check maps!"
572 << " " << board_id << " " << fEventTree->GetLeaf("tofpetId")->GetValue(n)
573 << " " << fEventTree->GetLeaf("tofpetChannel")->GetValue(n)
574 << " " << fEventTree->GetLeaf("tac")->GetValue(n)
575 << " " << fEventTree->GetLeaf("tCoarse")->GetValue(n)
576 << " " << fEventTree->GetLeaf("tFine")->GetValue(n)
577 << " " << fEventTree->GetLeaf("vCoarse")->GetValue(n)
578 << " " << fEventTree->GetLeaf("vFine")->GetValue(n);
579 }
580
581 t1 = high_resolution_clock::now();
582 if ( Chi2ndof > chi2Max )
583 {
584 if (QDC>1E20) QDC = 997.; // checking for inf
585 if (QDC != QDC) QDC = 998.; // checking for nan
586 if (QDC>0) QDC = -QDC;
587 mask = true;
588 }
589 else if (satur > saturationLimit || QDC>1E20 || QDC != QDC)
590 {
591 if (QDC>1E20) QDC = 987.; // checking for inf
592 LOG (info) << "inf " << board_id << " " << fEventTree->GetLeaf("tofpetId")->GetValue(n)
593 << " " << fEventTree->GetLeaf("tofpetChannel")->GetValue(n)
594 << " " << fEventTree->GetLeaf("tac")->GetValue(n)
595 << " " << fEventTree->GetLeaf("vCoarse")->GetValue(n)
596 << " " << fEventTree->GetLeaf("vFine")->GetValue(n)
597 << " " << TDC-fEventTree->GetLeaf("tCoarse")->GetValue(n)
598 << " " << eventNumber << " " << Chi2ndof;
599 if (QDC != QDC) QDC = 988.; // checking for nan
600 LOG (info) << "nan " << board_id << " " << fEventTree->GetLeaf("tofpetId")->GetValue(n)
601 << " " << fEventTree->GetLeaf("tofpetChannel")->GetValue(n)
602 << " " << fEventTree->GetLeaf("tac")->GetValue(n)
603 << " " << fEventTree->GetLeaf("vCoarse")->GetValue(n)
604 << " " << fEventTree->GetLeaf("vFine")->GetValue(n)
605 << " " << TDC-fEventTree->GetLeaf("tCoarse")->GetValue(n)
606 << " " << TDC << " " << fEventTree->GetLeaf("tCoarse")->GetValue(n)
607 << " " << eventNumber << " " << Chi2ndof;
608 A = int(min(QDC,double(1000.)));
609 B = min(satur,double(999.))/1000.;
610 QDC = A+B;
611 mask = true;
612 }
613 else if ( Chi2ndof > chi2Max )
614 {
615 if (QDC>0) QDC = -QDC;
616 mask = true;
617 }
618 LOG (info) << "calibrated: tdc = " << TDC << ", qdc = " << QDC;//TDC clock cycle = 160 MHz or 6.25ns
619 t4 = high_resolution_clock::now();
620 // Set the unit of the execution time measurement to ns
621 counters["qdc"]+= duration_cast<nanoseconds>(t1 - t0).count();
622 counters["make"]+= duration_cast<nanoseconds>(t4-t0).count();
623
624 // MuFilter encoding
625 if (!scifi)
626 {
627 system = MufiSystem[board_id][tofpet_id];
628 key = (tofpet_id%2)*1000 + tofpet_channel;
629 tmp = boardMapsMu["MuFilter"][board_name][slots[tofpet_id]];
630 if (debug || !(tmp.find("not") == string::npos))
631 {
632 LOG (info) << system << " " << key << " " << board_name << " " << tofpet_id
633 << " " << tofpet_id%2 << " " << tofpet_channel;
634 }
635 sipmChannel = 99;
636 if ( TofpetMap[system].count(key) == 0)
637 {
638 cout << "key " << key << " does not exist. " << endl
639 << board_name << " Tofpet id " << tofpet_id
640 << " System " << system << " has Tofpet map elements: {";
641 for ( auto it : TofpetMap[system])
642 {
643 cout << it.first << " : " << it.second << ", ";
644 }
645 cout << "}\n";
646 }
647 else sipmChannel = TofpetMap[system][key]-1;
648
649 nSiPMs = abs(offMap[tmp][1]);
650 nSides = abs(offMap[tmp][2]);
651 direction = int(offMap[tmp][1]/nSiPMs);
652 detID = offMap[tmp][0] + direction*int(sipmChannel/nSiPMs);
653 sipm_number = sipmChannel%(nSiPMs);
654 if ( tmp.find("Right") != string::npos ) sipm_number+= nSiPMs;
655 if (digiMuFilterStore.count(detID)==0)
656 {
657 digiMuFilterStore[detID] = new MuFilterHit(detID,nSiPMs,nSides);
658 }
659 test = digiMuFilterStore[detID]->GetSignal(sipm_number);
660 digiMuFilterStore[detID]->SetDigi(QDC,TDC,sipm_number);
661 digiMuFilterStore[detID]->SetDaqID(sipm_number,n, board_id, tofpet_id, tofpet_channel);
662 if (mask) digiMuFilterStore[detID]->SetMasked(sipm_number);
663
664 LOG (info) << "create mu hit: " << detID << " " << tmp << " " << system
665 << " " << tofpet_id << " " << offMap[tmp][0] << " " << offMap[tmp][1]
666 << " " << offMap[tmp][2] << " " << sipmChannel << " " << nSiPMs
667 << " " << nSides << " " << test << endl
668 << detID << " " << sipm_number << " " << QDC << " " << TDC;
669
670 if (test>0 || detID%1000>200 || sipm_number>15)
671 {
672 cout << "what goes wrong? " << detID << " SiPM " << sipm_number << " system " << system
673 << " key " << key << " board " << board_name << " tofperID " << tofpet_id
674 << " tofperChannel " << tofpet_channel << " test " << test << endl;
675 }
676 t5 = high_resolution_clock::now();
677 counters["createMufi"]+= duration_cast<nanoseconds>(t5 - t4).count();
678 } // end MuFilter encoding
679
680 else // now Scifi encoding
681 {
682 chan = channel_func(tofpet_id, tofpet_channel, mat);
683 orientation = 1;
684 if (station[2]=='Y') orientation = 0;
685 sipmLocal = (chan - mat*512);
686 sipmID = 1000000*int(station[1]-'0') + 100000*orientation + 10000*mat
687 + 1000*(int(sipmLocal/128)) + chan%128;
688 if (digiSciFiStore.count(sipmID)==0)
689 {
690 digiSciFiStore[sipmID] = new sndScifiHit(sipmID);
691 }
692 digiSciFiStore[sipmID]->SetDigi(QDC,TDC);
693 digiSciFiStore[sipmID]->SetDaqID(0,n,board_id, tofpet_id, tofpet_channel);
694 if (mask) digiSciFiStore[sipmID]->setInvalid();
695 LOG (info) << "create scifi hit: tdc = " << board_name << " " << sipmID
696 << " " << QDC << " " << TDC <<endl
697 << "tofpet:" << " " << tofpet_id << " " << tofpet_channel << " " << mat
698 << " " << chan << endl
699 << station[1] << " " << station[2] << " " << mat << " " << chan
700 << " " << int(chan/128)%4 << " " << chan%128;
701 t5 = high_resolution_clock::now();
702 counters["createScifi"]+= duration_cast<nanoseconds>(t5 - t4).count();
703 } // end Scifi encoding
704 } // end loop over hits in event
705
706 counters["N"]+= 1;
707 t6 = high_resolution_clock::now();
708 for (auto it_sipmID : digiSciFiStore)
709 {
710 (*fDigiSciFi)[indexSciFi]=digiSciFiStore[it_sipmID.first];
711 indexSciFi+= 1;
712 }
713 for (auto it_detID : digiMuFilterStore)
714 {
715 (*fDigiMuFilter)[indexMuFilter]=digiMuFilterStore[it_detID.first];
716 indexMuFilter+= 1;
717 }
718 counters["storage"]+= duration_cast<nanoseconds>(high_resolution_clock::now() - t6).count();
719 counters["event"]+= duration_cast<nanoseconds>(high_resolution_clock::now() - tE).count();
720 //timer.Stop();
721 //cout<<timer.RealTime()<<endl;
722
723 LOG (info) << fnStart+1 << " events processed out of "
724 << fEventTree->GetEntries() << " number of events in file.";
725 /*
726 cout << "Monitor computing time:" << endl;
727 for (auto it: counters)
728 {
729 if( it.first=="N" )
730 {
731 cout << "Processed " << it.first<< " = " << it.second << " events." << endl;
732 }
733 else
734 {
735 // Print execution time in seconds
736 cout << "stage " << it.first<< " took " << it.second*1e-9 << " [s]" << endl;
737 }
738 }
739 */
740}
741
742/* https://gitlab.cern.ch/snd-scifi/software/-/wikis/Raw-data-format
743 tac: 0-3, identifies the Time-To-Analogue converter used for this hit (each channel has four of them and they require separate calibration).
744 t_coarse: Coarse timestamp of the hit, running on a 4 times the LHC clock
745 t_fine: 0-1023, fine timestamp of the hit. It contains the raw value of the charge digitized by the TDC and requires calibration.
746 v_coarse: 0-1023, QDC mode: it represents the number of clock cycles the charge integration lasted.
747 v_fine = 0-1023, QDC mode: represents the charge measured. Requires calibration.
748*/
750void ConvRawData::StartTimeofRun(string Path)
751{
752 nlohmann::json j;
753 // reading file with xrootd protocol
754 File file;
755 XRootDStatus status;
756 StatInfo *info;
757 uint64_t offset = 0;
758 uint32_t size;
759 uint32_t bytesRead = 0;
760 if (local)
761 {
762 ifstream jsonfile(Form("%s/run_timestamps.json", Path.c_str()));
763 jsonfile >> j;
764 }
765 else
766 {
767 status = file.Open(Form("%s/run_timestamps.json", Path.c_str()), OpenFlags::Read);
768 if( !status.IsOK() )
769 {
770 LOG (error) << "Error opening file";
771 exit(0);
772 }
773 file.Stat(false, info);
774 size = info->GetSize();
775 char *buff = new char[size];
776 status = file.Read(offset, size, buff, bytesRead);
777 vector<char> vec;
778 for (size_t i = 0; i < size; i++){vec.push_back(buff[i]);}
779 j = json::parse(vec);
780 status = file.Close();
781 delete info;
782 delete [] buff;
783 }
784 // Read the start-time string
785 struct tm tm;
786 for (auto& el : j.items())
787 {
788 if (el.key() == "start_time")
789 {
790 string timeStr = el.value().get<string>();
791 strptime(timeStr.c_str(), "%Y-%m-%dT%H:%M:%SZ", &tm);
792 }
793 }
794 runStartUTC = timegm(&tm);
795}
797void ConvRawData::DetMapping(string Path)
798{
799 nlohmann::json j;
800 // reading file with xrootd protocol
801 File file;
802 XRootDStatus status;
803 StatInfo *info;
804 uint64_t offset = 0;
805 uint32_t size;
806 uint32_t bytesRead = 0;
807 // Call boardMappingParser
808 if (local)
809 {
810 ifstream jsonfile(Form("%s/board_mapping.json", Path.c_str()));
811 jsonfile >> j;
812 }
813 else
814 {
815 status = file.Open(Form("%s/board_mapping.json", Path.c_str()), OpenFlags::Read);
816 if( !status.IsOK() )
817 {
818 LOG (error) << "Error opening file";
819 exit(0);
820 }
821 file.Stat(false, info);
822 size = info->GetSize();
823 char *buff = new char[size];
824 status = file.Read(offset, size, buff, bytesRead);
825 vector<char> vec;
826 for (size_t i = 0; i < size; i++){vec.push_back(buff[i]);}
827 j = json::parse(vec);
828 status = file.Close();
829 delete info;
830 delete [] buff;
831 }
832 // Pass the read string to getBoardMapping()
833 tie(boardMaps, boardMapsMu) = getBoardMapping(j);
834
835 int board_id_mu{}, s{};
836 string tmp;
837 for ( auto it : boardMapsMu["MuFilter"] )
838 {
839 board_id_mu = stoi(it.first.substr(it.first.find("_") + 1));
840 for ( auto x : boardMapsMu["MuFilter"][it.first] )
841 {
842 for ( auto slot : slots )
843 {
844 if ( Path.find("testbeam_23") == string::npos ) s = 0;
845 else s = 3;
846 tmp = x.second.substr(0, x.second.find("_"));
847 if ( tmp =="US" ) s = 1;
848 if ( tmp=="DS" )
849 {
850 s = 2;
851 if ( Path.find("testbeam_24") != string::npos && x.second.substr(x.second.find("_")+1,1)==2) s = 3;
852 }
853 if ( slots[slot.first] == x.first )
854 {
855 MufiSystem[board_id_mu][slot.first] = s;
856 boardMaps["MuFilter"][it.first][slot.second] = slot.first;
857 }
858 }
859 }
860 }
861
862 //string o;
863 for (int i = 1 ; i < 6; i++ )
864 {
865 if (i<3)
866 {
867 offMap[Form("Veto_%iLeft",i)] = {10000 + (i-1)*1000+ 0, 8, 2};//first channel, nSiPMs, nSides
868 offMap[Form("Veto_%iRight",i)] = {10000 + (i-1)*1000+ 0, 8, 2};
869 }
870 if (i==3) offMap[Form("Veto_%iVert",i)] = {10000 + (i-1)*1000+ 6, -8, 1};// 3rd Veto plane
871 if (i<4)
872 {
873 // DS
874 offMap[Form("DS_%iLeft",i)] = {30000 + (i-1)*1000+ 59, -1, 2};// direction not known
875 offMap[Form("DS_%iRight",i)] = {30000 + (i-1)*1000+ 59, -1, 2};// direction not known
876 }
877 if (i<5)
878 {
879 offMap[Form("DS_%iVert",i)] = {30000 + (i-1)*1000+ 119, -1, 1};// direction not known
880 }
881 offMap[Form("US_%iLeft",i)] = {20000 + (i-1)*1000+ 9, -8, 2}; // from top to bottom
882 offMap[Form("US_%iRight",i)] = {20000 + (i-1)*1000+ 9, -8, 2};
883 }
884}
885/* Define calibration functions **/
886double ConvRawData::qdc_calibration( int board_id, int tofpet_id, int channel,
887 int tac, uint16_t v_coarse, uint16_t v_fine, uint16_t tf )
888 {
889 double GQDC = 1.0; // or 3.6
890 map<string, double> par= X_qdc[{board_id,tofpet_id,channel,tac}];
891 double x = v_coarse - tf;
892 double fqdc = -par["c"]*log(1+exp( par["a"]*pow((x-par["e"]),2)
893 -par["b"]*(x-par["e"]) )) + par["d"];
894 double value = (v_fine-fqdc)/GQDC;
895 return value;
896}
897double ConvRawData::qdc_chi2( int board_id, int tofpet_id, int channel, int tac, int TDC=0 )
898{
899 map<string, double> par = X_qdc[{board_id,tofpet_id,channel,tac}];
900 map<string, double> parT = X_tdc[{board_id,tofpet_id,channel,tac, TDC}];
901 return max(par["chi2Ndof"],parT["chi2Ndof"]);
902}
903double ConvRawData::qdc_sat( int board_id, int tofpet_id, int channel, int tac, uint16_t v_fine )
904{
905 map<string, double> par = X_qdc[{board_id,tofpet_id,channel,tac}];
906 return v_fine/par["d"];
907}
908double ConvRawData::time_calibration( int board_id, int tofpet_id, int channel,
909 int tac, int64_t t_coarse, uint16_t t_fine, int TDC=0 )
910{
911 map<string, double> parT = X_tdc[{board_id,tofpet_id,channel,tac, TDC}];
912 double x = t_fine;
913 double ftdc = (-parT["b"]-sqrt(pow(parT["b"],2)
914 -4*parT["a"]*(parT["c"]-x)))/(2*parT["a"]);
915 double timestamp = t_coarse+ftdc;
916 return timestamp;
917}
918tuple<double, double, double, double> ConvRawData::comb_calibration( int board_id, int tofpet_id, int channel,int tac, uint16_t v_coarse, uint16_t v_fine,int64_t t_coarse, uint16_t t_fine, double GQDC = 1.0, int TDC=0 )// max gain QDC = 3.6
919{
920 map<string, double> par= X_qdc[{board_id,tofpet_id,channel,tac}];
921 map<string, double> parT = X_tdc[{board_id,tofpet_id,channel,tac,TDC}];
922 double x = t_fine;
923 double ftdc = (-parT["b"]-sqrt(pow(parT["b"],2)
924 -4*parT["a"]*(parT["c"]-x)))/(2*parT["a"]); // Ettore 28/01/2022 +par['d']
925 double timestamp = t_coarse + ftdc;
926 double tf = timestamp - t_coarse;
927 x = v_coarse - tf;
928 double fqdc = - par["c"]*log(1+exp( par["a"]*pow((x-par["e"]),2)-par["b"]*(x-par["e"]) ))
929 + par["d"];
930 double value = (v_fine-fqdc)/GQDC;
931 return make_tuple(timestamp,value,max(par["chi2Ndof"],parT["chi2Ndof"]),v_fine/par["d"]);
932}
933map<double, pair<double, double> > ConvRawData::calibrationReport()
934{
935 TH1F* h = new TH1F("chi2","chi2", 1000, 0., 10000);
936 map<double, pair<double, double> > report{};
937 int TDC = 0;
938 double chi2{}, chi2T{}, key{};
939 // loop over entries of X_qdc and get map keys
940 for (auto it : X_qdc)
941 {
942 map<string, double> par= X_qdc[{it.first[0], it.first[1], it.first[2], it.first[3]}];
943 if (par["chi2Ndof"]) chi2 = par["chi2Ndof"];
944 else chi2 = -1;
945 map<string, double> parT = X_tdc[{it.first[0], it.first[1], it.first[2], it.first[3], TDC}];
946 if (parT["chi2Ndof"]) chi2T = parT["chi2Ndof"];
947 else chi2T = -1;
948 key = it.first[3] +10*it.first[2] + it.first[1]*10*100 + it.first[0]*10*100*100;
949 if (report.count(key)==0) report[key] = make_pair(chi2,chi2T);
950 }
951 for (auto it : report)
952 {
953 h->Fill(report[it.first].first);
954 h->Fill(report[it.first].second);
955 }
956 h->Draw();
957 return report;
958}
960int ConvRawData::channel_func( int tofpet_id, int tofpet_channel, int position)
961{
962 return (64*tofpet_id + 63 - tofpet_channel + 512*position);// 512 channels per mat, 1536 channels per plane. One channel covers all 6 layers of fibres.
963}
965void ConvRawData::read_csv(string Path)
966{
967 ifstream infile; // ifstream is the stream class to only read! from files.
968 stringstream X;
969 // reading file with xrootd protocol
970 File file;
971 XRootDStatus status;
972 StatInfo *info;
973 uint64_t offset = 0;
974 uint32_t size;
975 uint32_t bytesRead = 0;
976 string line, element;
977 double chi2_Ndof{};
978 // counter of number of elements
979 vector<int> key_vector{};
980
981 // Always read the SiPM mapping
982 int SiPM{};
983 vector<int> data_vector{};
984 map<string, map<int, vector<int>> > SiPMmap{};
985 vector<int> row{};
986 map<string, int> key { {"BM",3}, {"DS",2}, {"US",1}, {"Veto",0} };
987 struct stat buffer;
988 TString sndRoot = gSystem->Getenv("SNDSW_ROOT");
989 string sndswPath = sndRoot.Data();
990 string path_SiPMmap = Form("%s/geometry", sndswPath.c_str());
991 if (stat(path_SiPMmap.c_str(), &buffer) != 0)
992 {
993 LOG (error) << "Path "<< path_SiPMmap.c_str() << " does not exist!";
994 exit (0);
995 }
996 for (auto sys : key)
997 {
998 infile.open(Form("%s/%s_SiPM_mapping.csv", path_SiPMmap.c_str(), sys.first.c_str()));
999 X << infile.rdbuf();
1000 infile.close();
1001
1002 // Skip 1st line
1003 getline(X,line);
1004 LOG (info) << "In " << sys.first << " SiPM map file: " << line;
1005 // Read all other lines
1006 while (getline(X,line))
1007 {
1008 data_vector.clear();
1009 stringstream items(line);
1010 // first element in line is SiPM
1011 getline(items, element, ',');
1012 SiPM = stoi(element);
1013 // only taking first few elements of line
1014 while (data_vector.size()<4 && getline(items, element, ','))
1015 {
1016 data_vector.push_back(stoi(element));
1017 }
1018 SiPMmap[sys.first][SiPM] = data_vector;
1019 if (X.peek() == EOF) break;
1020 }
1021 X.str(string()); X.clear(); line.clear();
1022 size = 0; offset = 0; bytesRead = 0;
1023 for (auto channel : SiPMmap[sys.first])
1024 {
1025 row = channel.second;
1026 TofpetMap[sys.second][row.at(2)*1000+row.at(3)] = channel.first;
1027 }
1028 } // end filling SiPMmap and TofpetMap
1029
1030 // Get QDC/TDC calibration data if flag is set
1031 if (!makeCalibration) return;
1032
1033 // Get QDC calibration data
1034 if (local)
1035 {
1036 infile.open(Form("%s/qdc_cal.csv", Path.c_str()));
1037 X << infile.rdbuf();
1038 infile.close();
1039 }
1040 else
1041 {
1042 status = file.Open(Form("%s/qdc_cal.csv", Path.c_str()), OpenFlags::Read);
1043 if( !status.IsOK() )
1044 {
1045 LOG (error) << "Error opening file";
1046 exit(0);
1047 }
1048 file.Stat(false, info);
1049 size = info->GetSize();
1050 char *buff = new char[size];
1051 status = file.Read(offset, size, buff, bytesRead);
1052 X << buff;
1053 status = file.Close();
1054 delete info;
1055 delete [] buff;
1056 }
1057 LOG (info) << "Read_csv "<<Path.c_str();
1058 // Skip 1st line
1059 getline(X, line);
1060 LOG (info) << "In QDC cal file: " << line;
1061 vector<double> qdcData{};
1062 // Read all other lines
1063 while (getline(X, line))
1064 {
1065 qdcData.clear();
1066 stringstream items(line);
1067 while (getline(items, element, ','))
1068 {
1069 if(iscntrl(element[0])) break;
1070 qdcData.push_back(stof(element));
1071 }
1072 if (qdcData.size()<10) continue;
1073 key_vector.clear();
1074 key_vector = {int(qdcData[0]), int(qdcData[1]), int(qdcData[2]), int(qdcData[3])};
1075 chi2_Ndof = (qdcData[9] < 2) ? 999999. : qdcData[7]/qdcData[9];
1076 X_qdc[key_vector] = { {"a",qdcData[4]}, {"b",qdcData[5]}, {"c",qdcData[6]},
1077 {"d",qdcData[8]}, {"e",qdcData[10]}, {"chi2Ndof",chi2_Ndof} };
1078 if (X.peek() == EOF) break;
1079 }
1080 X.str(string()); X.clear(); line.clear();
1081 size = 0; offset = 0; bytesRead = 0;
1082 // Get TDC calibration data
1083 if (local)
1084 {
1085 infile.open(Form("%s/tdc_cal.csv", Path.c_str()));
1086 X << infile.rdbuf();
1087 infile.close();
1088 }
1089 else
1090 {
1091 status = file.Open(Form("%s/tdc_cal.csv", Path.c_str()), OpenFlags::Read);
1092 if( !status.IsOK() )
1093 {
1094 LOG (error) << "Error opening file";
1095 exit(0);
1096 }
1097 file.Stat(false, info);
1098 size = info->GetSize();
1099 char *buff = new char[size];
1100 status = file.Read(offset, size, buff, bytesRead);
1101 X << buff;
1102 status = file.Close();
1103 delete info;
1104 delete [] buff;
1105 }
1106 // Skip 1st line
1107 getline(X, line);
1108 LOG (info) << "In TDC cal file: " << line;
1109 vector<double> tdcData{};
1110 // Read all other lines
1111 while (getline(X, line))
1112 {
1113 tdcData.clear();
1114 stringstream items(line);
1115 if (line.length()<5){continue;}
1116 while (getline(items, element, ','))
1117 {
1118 if(iscntrl(element[0])) break;
1119 tdcData.push_back(stof(element));
1120 }
1121 if (tdcData.size()<9) continue;
1122 key_vector.clear();
1123 key_vector = {int(tdcData[0]), int(tdcData[1]), int(tdcData[2]), int(tdcData[3]), int(tdcData[4])};
1124 chi2_Ndof = (tdcData[10] < 2) ? 999999. : tdcData[8]/tdcData[10];
1125 X_tdc[key_vector] = { {"a",tdcData[5]}, {"b",tdcData[6]}, {"c",tdcData[7]},
1126 {"d",tdcData[9]}, {"chi2Ndof",chi2_Ndof} };
1127 if (X.peek() == EOF) break;
1128 }
1129 X.str(string()); X.clear(); line.clear();
1130 size = 0; offset = 0; bytesRead = 0;
1131}
1132void ConvRawData::debugMapping(string brd, int tofpetID, int tofpetChannel)
1133{
1134 int Key{}, SiPMChannel{}, n_SiPMs{}, n_Sides{}, Direction{}, DetID{}, SiPM_number{};
1135 string Tmp;
1136 int brdID{}, System{};
1137 Key = (tofpetID%2)*1000 + tofpetChannel;
1138 Tmp = boardMapsMu["MuFilter"][brd][slots[tofpetID]];
1139 brdID = stoi(brd.substr(brd.find("_")+1));
1140 System = MufiSystem[brdID][tofpetID];
1141 SiPMChannel = TofpetMap[System][Key]-1;
1142 n_SiPMs = abs(offMap[Tmp][1]);
1143 n_Sides = abs(offMap[Tmp][2]);
1144 Direction = int(offMap[Tmp][1]/n_SiPMs);
1145 DetID = offMap[Tmp][0] + Direction*int(SiPMChannel/n_SiPMs);
1146 SiPM_number = SiPMChannel%(n_SiPMs);
1147 cout << "SiPM channel " << SiPMChannel << " nSiPM " << n_SiPMs << " nSides " << n_Sides
1148 << " detID " << DetID << " SiPM number " << SiPM_number << endl;
1149}
1150
1151ClassImp(ConvRawData);
ClassImp
ClassImp(ConvRawData)
Direction
Direction
Definition RPCUnpack.cxx:23
boardMaps
map< string, map< string, map< string, int > > > boardMaps
Definition boardMappingParser.cxx:14
getBoardMapping
tuple< map< string, map< string, map< string, int > > >, map< string, map< string, map< string, string > > > > getBoardMapping(json j)
Definition boardMappingParser.cxx:65
boardMapsMu
map< string, map< string, map< string, string > > > boardMapsMu
Definition boardMappingParser.cxx:15
ConvRawData::qdc_chi2
double qdc_chi2(int board_id, int tofpet_id, int channel, int tac, int TDC)
Definition ConvRawData.cxx:897
ConvRawData::fEventHeader
FairEventHeader * fEventHeader
Definition ConvRawData.h:104
ConvRawData::debugMapping
void debugMapping(std::string board, int tofpetID, int tofpetChannel)
Definition ConvRawData.cxx:1132
ConvRawData::qdc_calibration
double qdc_calibration(int board_id, int tofpet_id, int channel, int tac, uint16_t v_coarse, uint16_t v_fine, uint16_t tf)
Definition ConvRawData.cxx:886
ConvRawData::time_calibration
double time_calibration(int board_id, int tofpet_id, int channel, int tac, int64_t t_coarse, uint16_t t_fine, int TDC)
Definition ConvRawData.cxx:908
ConvRawData::counters
std::map< std::string, double > counters
Definition ConvRawData.h:81
ConvRawData::Exec
virtual void Exec(Option_t *opt)
Definition ConvRawData.cxx:153
ConvRawData::fDigiMuFilter
TClonesArray * fDigiMuFilter
Definition ConvRawData.h:106
ConvRawData::offMap
std::map< std::string, std::vector< int > > offMap
Definition ConvRawData.h:79
ConvRawData::DetMapping
void DetMapping(std::string path)
Definition ConvRawData.cxx:797
ConvRawData::Init
virtual InitStatus Init()
Definition ConvRawData.cxx:60
ConvRawData::ScifiDet
Scifi * ScifiDet
Definition ConvRawData.h:84
ConvRawData::fpathCalib
std::string fpathCalib
Definition ConvRawData.h:95
ConvRawData::saturationLimit
double saturationLimit
Definition ConvRawData.h:97
ConvRawData::StartTimeofRun
void StartTimeofRun(std::string path)
Definition ConvRawData.cxx:750
ConvRawData::comb_calibration
std::tuple< double, double, double, double > comb_calibration(int board_id, int tofpet_id, int channel, int tac, uint16_t v_coarse, uint16_t v_fine, int64_t t_coarse, uint16_t t_fine, double GQDC, int TDC)
Definition ConvRawData.cxx:918
ConvRawData::digiMuFilterStore
std::map< int, MuFilterHit * > digiMuFilterStore
Definition ConvRawData.h:69
ConvRawData::slots
std::map< int, std::string > slots
Definition ConvRawData.h:75
ConvRawData::UpdateInput
void UpdateInput(int n)
Definition ConvRawData.cxx:166
ConvRawData::chi2Max
double chi2Max
Definition ConvRawData.h:97
ConvRawData::fDigiSciFi
TClonesArray * fDigiSciFi
Definition ConvRawData.h:105
ConvRawData::MufiSystem
std::map< int, std::map< int, int > > MufiSystem
Definition ConvRawData.h:74
ConvRawData::qdc_sat
double qdc_sat(int board_id, int tofpet_id, int channel, int tac, uint16_t v_fine)
Definition ConvRawData.cxx:903
ConvRawData::X_qdc
std::map< std::vector< int >, std::map< std::string, double > > X_qdc
Definition ConvRawData.h:71
ConvRawData::fpathJSON
std::string fpathJSON
Definition ConvRawData.h:95
ConvRawData::fSNDLHCEventHeader
SNDLHCEventHeader * fSNDLHCEventHeader
Definition ConvRawData.h:103
ConvRawData::fEventTree
TTree * fEventTree
Definition ConvRawData.h:87
ConvRawData::runStartUTC
double runStartUTC
Definition ConvRawData.h:98
ConvRawData::boards
std::map< std::string, TTree * > boards
Definition ConvRawData.h:89
ConvRawData::X_tdc
std::map< std::vector< int >, std::map< std::string, double > > X_tdc
Definition ConvRawData.h:72
ConvRawData::digiSciFiStore
std::map< int, sndScifiHit * > digiSciFiStore
Definition ConvRawData.h:70
ConvRawData::channel_func
int channel_func(int tofpet_id, int tofpet_channel, int position)
Definition ConvRawData.cxx:960
ConvRawData::TofpetMap
std::map< int, std::map< int, int > > TofpetMap
Definition ConvRawData.h:77
ConvRawData::read_csv
void read_csv(std::string path)
Definition ConvRawData.cxx:965
ConvRawData::makeCalibration
int makeCalibration
Definition ConvRawData.h:94
ConvRawData::calibrationReport
std::map< double, std::pair< double, double > > calibrationReport()
Definition ConvRawData.cxx:933
SNDLHCEventHeader::SetBunchType
void SetBunchType(int16_t bunchType)
Definition SNDLHCEventHeader.h:35
SNDLHCEventHeader::SetEventTime
void SetEventTime(int64_t time)
Definition SNDLHCEventHeader.h:30
SNDLHCEventHeader::SetRunId
void SetRunId(uint64_t runid)
Definition SNDLHCEventHeader.h:29
SNDLHCEventHeader::SetFlags
void SetFlags(uint64_t flags)
Definition SNDLHCEventHeader.cxx:60
SNDLHCEventHeader::SetUTCtimestamp
void SetUTCtimestamp(int64_t UTCtstamp)
Definition SNDLHCEventHeader.h:33
SNDLHCEventHeader::SetEventNumber
void SetEventNumber(int id)
Definition SNDLHCEventHeader.h:32
Scifi
Definition Scifi.h:20