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 // mini DTs are labelled DS 5V, board is 32
337 if (tmp == "DS_5Vert") continue;
338 if (debug || !(tmp.find("not") == string::npos))
339 {
340 LOG (info) << system << " " << key << " " << board.first << " " << tofpet_id
341 << " " << tofpet_id%2 << " " << tofpet_channel;
342 }
343 sipmChannel = 99;
344 if ( TofpetMap[system].count(key) == 0)
345 {
346 cout << "key " << key << " does not exist. " << endl
347 << board.first << " Tofpet id " << tofpet_id
348 << " System " << system << " has Tofpet map elements: {";
349 for ( auto it : TofpetMap[system])
350 {
351 cout << it.first << " : " << it.second << ", ";
352 }
353 cout << "}\n";
354 }
355 else sipmChannel = TofpetMap[system][key]-1;
356
357 nSiPMs = abs(offMap[tmp][1]);
358 nSides = abs(offMap[tmp][2]);
359 direction = int(offMap[tmp][1]/nSiPMs);
360 detID = offMap[tmp][0] + direction*int(sipmChannel/nSiPMs);
361 sipm_number = sipmChannel%(nSiPMs);
362 if ( tmp.find("Right") != string::npos ) sipm_number+= nSiPMs;
363 if (digiMuFilterStore.count(detID)==0)
364 {
365 digiMuFilterStore[detID] = new MuFilterHit(detID,nSiPMs,nSides);
366 }
367 test = digiMuFilterStore[detID]->GetSignal(sipm_number);
368 digiMuFilterStore[detID]->SetDigi(QDC,TDC,sipm_number);
369 digiMuFilterStore[detID]->SetDaqID(sipm_number,n, board_id, tofpet_id, tofpet_channel);
370 if (mask) digiMuFilterStore[detID]->SetMasked(sipm_number);
371
372 LOG (info) << "create mu hit: " << detID << " " << tmp << " " << system
373 << " " << tofpet_id << " " << offMap[tmp][0] << " " << offMap[tmp][1]
374 << " " << offMap[tmp][2] << " " << sipmChannel << " " << nSiPMs
375 << " " << nSides << " " << test << endl
376 << detID << " " << sipm_number << " " << QDC << " " << TDC;
377
378 if (test>0 || detID%1000>200 || sipm_number>15)
379 {
380 cout << "what goes wrong? " << detID << " SiPM " << sipm_number << " system " << system
381 << " key " << key << " board " << board.first << " tofperID " << tofpet_id
382 << " tofperChannel " << tofpet_channel << " test " << test << endl;
383 }
384 t5 = high_resolution_clock::now();
385 counters["createMufi"]+= duration_cast<nanoseconds>(t5 - t4).count();
386 } // end MuFilter encoding
387
388 else // now Scifi encoding
389 {
390 chan = channel_func(tofpet_id, tofpet_channel, mat);
391 orientation = 1;
392 if (station[2]=='Y') orientation = 0;
393 sipmLocal = (chan - mat*512);
394 sipmID = 1000000*int(station[1]-'0') + 100000*orientation + 10000*mat
395 + 1000*(int(sipmLocal/128)) + chan%128;
396 if (digiSciFiStore.count(sipmID)==0)
397 {
398 digiSciFiStore[sipmID] = new sndScifiHit(sipmID);
399 }
400 digiSciFiStore[sipmID]->SetDigi(QDC,TDC);
401 digiSciFiStore[sipmID]->SetDaqID(0,n, board_id, tofpet_id, tofpet_channel);
402 if (mask) digiSciFiStore[sipmID]->setInvalid();
403 LOG (info) << "create scifi hit: tdc = " << board.first << " " << sipmID
404 << " " << QDC << " " << TDC <<endl
405 << "tofpet:" << " " << tofpet_id << " " << tofpet_channel << " " << mat
406 << " " << chan << endl
407 << station[1] << " " << station[2] << " " << mat << " " << chan
408 << " " << int(chan/128)%4 << " " << chan%128;
409 t5 = high_resolution_clock::now();
410 counters["createScifi"]+= duration_cast<nanoseconds>(t5 - t4).count();
411 } // end Scifi encoding
412 } // end loop over hits in event
413 } // end loop over boards (TTrees in data file)
414
415 counters["N"]+= 1;
416 t6 = high_resolution_clock::now();
417 for (auto it_sipmID : digiSciFiStore)
418 {
419 (*fDigiSciFi)[indexSciFi]=digiSciFiStore[it_sipmID.first];
420 indexSciFi+= 1;
421 }
422 for (auto it_detID : digiMuFilterStore)
423 {
424 (*fDigiMuFilter)[indexMuFilter]=digiMuFilterStore[it_detID.first];
425 indexMuFilter+= 1;
426 }
427 counters["storage"]+= duration_cast<nanoseconds>(high_resolution_clock::now() - t6).count();
428 counters["event"]+= duration_cast<nanoseconds>(high_resolution_clock::now() - tE).count();
429 //timer.Stop();
430 //cout<<timer.RealTime()<<endl;
431
432 LOG (info) << fnStart+1 << " events processed out of "
433 << fEventTree->GetEntries() << " number of events in file.";
434 /*
435 cout << "Monitor computing time:" << endl;
436 for (auto it: counters)
437 {
438 if( it.first=="N" )
439 {
440 cout << "Processed " << it.first<< " = " << it.second << " events." << endl;
441 }
442 else
443 {
444 // Print execution time in seconds
445 cout << "stage " << it.first<< " took " << it.second*1e-9 << " [s]" << endl;
446 }
447 }
448 */
449}
450
451void ConvRawData::Process1()
452{
453 int indexSciFi{}, indexMuFilter{};
454 bool scifi = false, mask = false;
455 int board_id{};
456 string board_name;
457 int tofpet_id{}, tofpet_channel{}, tac{}, mat{};
458 string station;
459 double TDC{}, QDC{}, Chi2ndof{}, satur{};
460 int A{};
461 double B{};
462 high_resolution_clock::time_point tE{}, t0{}, t1{}, t4{}, t5{}, t6{}, tt{};
463 int system{}, key{}, sipmChannel{};
464 string tmp;
465 int nSiPMs{}, nSides{}, direction{}, detID{}, sipm_number{}, chan{}, orientation{}, sipmLocal{};
466 int sipmID{};
467 double test{};
468 int64_t eventTime{};
469 //TStopwatch timer;
470
471 tE = high_resolution_clock::now();
472 //timer.Start();
473 fEventTree->GetEvent(eventNumber);
474 if ( eventNumber%fheartBeat == 0 )
475 {
476 tt = high_resolution_clock::now();
477 time_t ttp = high_resolution_clock::to_time_t(tt);
478 LOG (info) << "run " << frunNumber << " event " << eventNumber
479 << " local time " << ctime(&ttp);
480 }
481
482 fSNDLHCEventHeader->SetFlags(fEventTree->GetLeaf("evtFlags")->GetValue());
483 fSNDLHCEventHeader->SetRunId(frunNumber);
484 eventTime = fEventTree->GetLeaf("evtTimestamp")->GetValue();
485 fSNDLHCEventHeader->SetEventTime(eventTime);
486 fSNDLHCEventHeader->SetUTCtimestamp(eventTime*6.23768*1e-9 + runStartUTC);
487 fSNDLHCEventHeader->SetEventNumber(fEventTree->GetLeaf("evtNumber")->GetValue());
488 // Fill filling scheme data into the event header
489 if (FSmap->GetEntries()>1){
490 // ion runs
491 if (fSNDLHCEventHeader->GetAccMode()== 12){
492 fSNDLHCEventHeader->SetBunchType(stoi(((TObjString*)FSmap->GetValue(Form("%d",
493 int((eventTime%(4*3564))/8+0.5))))->GetString().Data()));
494 }
495 // proton runs
496 else{
497 fSNDLHCEventHeader->SetBunchType(stoi(((TObjString*)FSmap->GetValue(Form("%d",
498 int((eventTime%(4*3564))/4))))->GetString().Data()));
499 }
500 }
501 else
502 fSNDLHCEventHeader->SetBunchType(stoi(((TObjString*)FSmap->GetValue("0"))->GetString().Data()));
503
504 LOG (info) << "evtNumber per run "
505 << fEventTree->GetLeaf("evtNumber")->GetValue()
506 << " evtNumber per partition: " << eventNumber
507 << " timestamp: " << eventTime;
508 // Delete pointer map elements
509 for (auto it : digiSciFiStore)
510 {
511 delete it.second;
512 }
513 digiSciFiStore.clear();
514 for (auto it : digiMuFilterStore)
515 {
516 delete it.second;
517 }
518 digiMuFilterStore.clear();
519
520 // Loop over hits per event!
521 for ( int n =0; n < fEventTree->GetLeaf("nHits")->GetValue(); n++ )
522 {
523 board_id = fEventTree->GetLeaf("boardId")->GetValue(n);
524 board_name = "board_"+to_string(board_id);
525 scifi = true;
526 if (boardMaps["Scifi"].count(board_name)!=0)
527 {
528 for (auto it : boardMaps["Scifi"][board_name])
529 {
530 station = it.first;
531 mat = it.second;
532 }
533 }
534 else if (boardMaps["MuFilter"].count(board_name)!=0) scifi = false;
535 else
536 {
537 LOG (error) << board_name << " not known. Serious error, stop!";
538 break;
539 }
540 mask = false;
541 LOG (info) << "In scifi? " << scifi
542 << " " << board_id << " " << fEventTree->GetLeaf("tofpetId")->GetValue(n)
543 << " " << fEventTree->GetLeaf("tofpetChannel")->GetValue(n)
544 << " " << fEventTree->GetLeaf("tac")->GetValue(n)
545 << " " << fEventTree->GetLeaf("tCoarse")->GetValue(n)
546 << " " << fEventTree->GetLeaf("tFine")->GetValue(n)
547 << " " << fEventTree->GetLeaf("vCoarse")->GetValue(n)
548 << " " << fEventTree->GetLeaf("vFine")->GetValue(n);
549 t0 = high_resolution_clock::now();
550 tofpet_id = fEventTree->GetLeaf("tofpetId")->GetValue(n);
551 tofpet_channel = fEventTree->GetLeaf("tofpetChannel")->GetValue(n);
552 tac = fEventTree->GetLeaf("tac")->GetValue(n);
553 /* Since run 39 calibration is done online and
554 calib data stored in root raw data file */
555 if (makeCalibration)
556 tie(TDC,QDC,Chi2ndof,satur) = comb_calibration(board_id, tofpet_id, tofpet_channel, tac,
557 fEventTree->GetLeaf("vCoarse")->GetValue(n),
558 fEventTree->GetLeaf("vFine")->GetValue(n),
559 fEventTree->GetLeaf("tCoarse")->GetValue(n),
560 fEventTree->GetLeaf("tFine")->GetValue(n),
561 1.0, 0);
562 else
563 {
564 TDC = fEventTree->GetLeaf("timestamp")->GetValue(n);
565 QDC = fEventTree->GetLeaf("value")->GetValue(n);
566 Chi2ndof = max(fEventTree->GetLeaf("timestampCalChi2")->GetValue(n)/fEventTree->GetLeaf("timestampCalDof")->GetValue(n),
567 fEventTree->GetLeaf("valueCalChi2")->GetValue(n)/fEventTree->GetLeaf("valueCalDof")->GetValue(n));
568 satur = fEventTree->GetLeaf("valueSaturation")->GetValue(n);
569 }
570
571 // Print a warning if TDC or QDC is nan.
572 if ( TDC != TDC || QDC!=QDC) {
573 LOG (error) << "NAN tdc/qdc detected! Check maps!"
574 << " " << board_id << " " << fEventTree->GetLeaf("tofpetId")->GetValue(n)
575 << " " << fEventTree->GetLeaf("tofpetChannel")->GetValue(n)
576 << " " << fEventTree->GetLeaf("tac")->GetValue(n)
577 << " " << fEventTree->GetLeaf("tCoarse")->GetValue(n)
578 << " " << fEventTree->GetLeaf("tFine")->GetValue(n)
579 << " " << fEventTree->GetLeaf("vCoarse")->GetValue(n)
580 << " " << fEventTree->GetLeaf("vFine")->GetValue(n);
581 }
582
583 t1 = high_resolution_clock::now();
584 if ( Chi2ndof > chi2Max )
585 {
586 if (QDC>1E20) QDC = 997.; // checking for inf
587 if (QDC != QDC) QDC = 998.; // checking for nan
588 if (QDC>0) QDC = -QDC;
589 mask = true;
590 }
591 else if (satur > saturationLimit || QDC>1E20 || QDC != QDC)
592 {
593 if (QDC>1E20) QDC = 987.; // checking for inf
594 LOG (info) << "inf " << board_id << " " << fEventTree->GetLeaf("tofpetId")->GetValue(n)
595 << " " << fEventTree->GetLeaf("tofpetChannel")->GetValue(n)
596 << " " << fEventTree->GetLeaf("tac")->GetValue(n)
597 << " " << fEventTree->GetLeaf("vCoarse")->GetValue(n)
598 << " " << fEventTree->GetLeaf("vFine")->GetValue(n)
599 << " " << TDC-fEventTree->GetLeaf("tCoarse")->GetValue(n)
600 << " " << eventNumber << " " << Chi2ndof;
601 if (QDC != QDC) QDC = 988.; // checking for nan
602 LOG (info) << "nan " << board_id << " " << fEventTree->GetLeaf("tofpetId")->GetValue(n)
603 << " " << fEventTree->GetLeaf("tofpetChannel")->GetValue(n)
604 << " " << fEventTree->GetLeaf("tac")->GetValue(n)
605 << " " << fEventTree->GetLeaf("vCoarse")->GetValue(n)
606 << " " << fEventTree->GetLeaf("vFine")->GetValue(n)
607 << " " << TDC-fEventTree->GetLeaf("tCoarse")->GetValue(n)
608 << " " << TDC << " " << fEventTree->GetLeaf("tCoarse")->GetValue(n)
609 << " " << eventNumber << " " << Chi2ndof;
610 A = int(min(QDC,double(1000.)));
611 B = min(satur,double(999.))/1000.;
612 QDC = A+B;
613 mask = true;
614 }
615 else if ( Chi2ndof > chi2Max )
616 {
617 if (QDC>0) QDC = -QDC;
618 mask = true;
619 }
620 LOG (info) << "calibrated: tdc = " << TDC << ", qdc = " << QDC;//TDC clock cycle = 160 MHz or 6.25ns
621 t4 = high_resolution_clock::now();
622 // Set the unit of the execution time measurement to ns
623 counters["qdc"]+= duration_cast<nanoseconds>(t1 - t0).count();
624 counters["make"]+= duration_cast<nanoseconds>(t4-t0).count();
625
626 // MuFilter encoding
627 if (!scifi)
628 {
629 system = MufiSystem[board_id][tofpet_id];
630 key = (tofpet_id%2)*1000 + tofpet_channel;
631 tmp = boardMapsMu["MuFilter"][board_name][slots[tofpet_id]];
632 // mini DTs are labelled DS 5V, board is 32
633 if (tmp == "DS_5Vert") continue;
634 if (debug || !(tmp.find("not") == string::npos))
635 {
636 LOG (info) << system << " " << key << " " << board_name << " " << tofpet_id
637 << " " << tofpet_id%2 << " " << tofpet_channel;
638 }
639 sipmChannel = 99;
640 if ( TofpetMap[system].count(key) == 0)
641 {
642 cout << "key " << key << " does not exist. " << endl
643 << board_name << " Tofpet id " << tofpet_id
644 << " System " << system << " has Tofpet map elements: {";
645 for ( auto it : TofpetMap[system])
646 {
647 cout << it.first << " : " << it.second << ", ";
648 }
649 cout << "}\n";
650 }
651 else sipmChannel = TofpetMap[system][key]-1;
652
653 nSiPMs = abs(offMap[tmp][1]);
654 nSides = abs(offMap[tmp][2]);
655 direction = int(offMap[tmp][1]/nSiPMs);
656 detID = offMap[tmp][0] + direction*int(sipmChannel/nSiPMs);
657 sipm_number = sipmChannel%(nSiPMs);
658 if ( tmp.find("Right") != string::npos ) sipm_number+= nSiPMs;
659 if (digiMuFilterStore.count(detID)==0)
660 {
661 digiMuFilterStore[detID] = new MuFilterHit(detID,nSiPMs,nSides);
662 }
663 test = digiMuFilterStore[detID]->GetSignal(sipm_number);
664 digiMuFilterStore[detID]->SetDigi(QDC,TDC,sipm_number);
665 digiMuFilterStore[detID]->SetDaqID(sipm_number,n, board_id, tofpet_id, tofpet_channel);
666 if (mask) digiMuFilterStore[detID]->SetMasked(sipm_number);
667
668 LOG (info) << "create mu hit: " << detID << " " << tmp << " " << system
669 << " " << tofpet_id << " " << offMap[tmp][0] << " " << offMap[tmp][1]
670 << " " << offMap[tmp][2] << " " << sipmChannel << " " << nSiPMs
671 << " " << nSides << " " << test << endl
672 << detID << " " << sipm_number << " " << QDC << " " << TDC;
673
674 if (test>0 || detID%1000>200 || sipm_number>15)
675 {
676 cout << "what goes wrong? " << detID << " SiPM " << sipm_number << " system " << system
677 << " key " << key << " board " << board_name << " tofperID " << tofpet_id
678 << " tofperChannel " << tofpet_channel << " test " << test << endl;
679 }
680 t5 = high_resolution_clock::now();
681 counters["createMufi"]+= duration_cast<nanoseconds>(t5 - t4).count();
682 } // end MuFilter encoding
683
684 else // now Scifi encoding
685 {
686 chan = channel_func(tofpet_id, tofpet_channel, mat);
687 orientation = 1;
688 if (station[2]=='Y') orientation = 0;
689 sipmLocal = (chan - mat*512);
690 sipmID = 1000000*int(station[1]-'0') + 100000*orientation + 10000*mat
691 + 1000*(int(sipmLocal/128)) + chan%128;
692 if (digiSciFiStore.count(sipmID)==0)
693 {
694 digiSciFiStore[sipmID] = new sndScifiHit(sipmID);
695 }
696 digiSciFiStore[sipmID]->SetDigi(QDC,TDC);
697 digiSciFiStore[sipmID]->SetDaqID(0,n,board_id, tofpet_id, tofpet_channel);
698 if (mask) digiSciFiStore[sipmID]->setInvalid();
699 LOG (info) << "create scifi hit: tdc = " << board_name << " " << sipmID
700 << " " << QDC << " " << TDC <<endl
701 << "tofpet:" << " " << tofpet_id << " " << tofpet_channel << " " << mat
702 << " " << chan << endl
703 << station[1] << " " << station[2] << " " << mat << " " << chan
704 << " " << int(chan/128)%4 << " " << chan%128;
705 t5 = high_resolution_clock::now();
706 counters["createScifi"]+= duration_cast<nanoseconds>(t5 - t4).count();
707 } // end Scifi encoding
708 } // end loop over hits in event
709
710 counters["N"]+= 1;
711 t6 = high_resolution_clock::now();
712 for (auto it_sipmID : digiSciFiStore)
713 {
714 (*fDigiSciFi)[indexSciFi]=digiSciFiStore[it_sipmID.first];
715 indexSciFi+= 1;
716 }
717 for (auto it_detID : digiMuFilterStore)
718 {
719 (*fDigiMuFilter)[indexMuFilter]=digiMuFilterStore[it_detID.first];
720 indexMuFilter+= 1;
721 }
722 counters["storage"]+= duration_cast<nanoseconds>(high_resolution_clock::now() - t6).count();
723 counters["event"]+= duration_cast<nanoseconds>(high_resolution_clock::now() - tE).count();
724 //timer.Stop();
725 //cout<<timer.RealTime()<<endl;
726
727 LOG (info) << fnStart+1 << " events processed out of "
728 << fEventTree->GetEntries() << " number of events in file.";
729 /*
730 cout << "Monitor computing time:" << endl;
731 for (auto it: counters)
732 {
733 if( it.first=="N" )
734 {
735 cout << "Processed " << it.first<< " = " << it.second << " events." << endl;
736 }
737 else
738 {
739 // Print execution time in seconds
740 cout << "stage " << it.first<< " took " << it.second*1e-9 << " [s]" << endl;
741 }
742 }
743 */
744}
745
746/* https://gitlab.cern.ch/snd-scifi/software/-/wikis/Raw-data-format
747 tac: 0-3, identifies the Time-To-Analogue converter used for this hit (each channel has four of them and they require separate calibration).
748 t_coarse: Coarse timestamp of the hit, running on a 4 times the LHC clock
749 t_fine: 0-1023, fine timestamp of the hit. It contains the raw value of the charge digitized by the TDC and requires calibration.
750 v_coarse: 0-1023, QDC mode: it represents the number of clock cycles the charge integration lasted.
751 v_fine = 0-1023, QDC mode: represents the charge measured. Requires calibration.
752*/
754void ConvRawData::StartTimeofRun(string Path)
755{
756 nlohmann::json j;
757 // reading file with xrootd protocol
758 File file;
759 XRootDStatus status;
760 StatInfo *info;
761 uint64_t offset = 0;
762 uint32_t size;
763 uint32_t bytesRead = 0;
764 if (local)
765 {
766 ifstream jsonfile(Form("%s/run_timestamps.json", Path.c_str()));
767 jsonfile >> j;
768 }
769 else
770 {
771 status = file.Open(Form("%s/run_timestamps.json", Path.c_str()), OpenFlags::Read);
772 if( !status.IsOK() )
773 {
774 LOG (error) << "Error opening file";
775 exit(0);
776 }
777 file.Stat(false, info);
778 size = info->GetSize();
779 char *buff = new char[size];
780 status = file.Read(offset, size, buff, bytesRead);
781 vector<char> vec;
782 for (size_t i = 0; i < size; i++){vec.push_back(buff[i]);}
783 j = json::parse(vec);
784 status = file.Close();
785 delete info;
786 delete [] buff;
787 }
788 // Read the start-time string
789 struct tm tm;
790 for (auto& el : j.items())
791 {
792 if (el.key() == "start_time")
793 {
794 string timeStr = el.value().get<string>();
795 strptime(timeStr.c_str(), "%Y-%m-%dT%H:%M:%SZ", &tm);
796 }
797 }
798 runStartUTC = timegm(&tm);
799}
801void ConvRawData::DetMapping(string Path)
802{
803 nlohmann::json j;
804 // reading file with xrootd protocol
805 File file;
806 XRootDStatus status;
807 StatInfo *info;
808 uint64_t offset = 0;
809 uint32_t size;
810 uint32_t bytesRead = 0;
811 // Call boardMappingParser
812 if (local)
813 {
814 ifstream jsonfile(Form("%s/board_mapping.json", Path.c_str()));
815 jsonfile >> j;
816 }
817 else
818 {
819 status = file.Open(Form("%s/board_mapping.json", Path.c_str()), OpenFlags::Read);
820 if( !status.IsOK() )
821 {
822 LOG (error) << "Error opening file";
823 exit(0);
824 }
825 file.Stat(false, info);
826 size = info->GetSize();
827 char *buff = new char[size];
828 status = file.Read(offset, size, buff, bytesRead);
829 vector<char> vec;
830 for (size_t i = 0; i < size; i++){vec.push_back(buff[i]);}
831 j = json::parse(vec);
832 status = file.Close();
833 delete info;
834 delete [] buff;
835 }
836 // Pass the read string to getBoardMapping()
837 tie(boardMaps, boardMapsMu) = getBoardMapping(j);
838
839 int board_id_mu{}, s{};
840 string tmp;
841 for ( auto it : boardMapsMu["MuFilter"] )
842 {
843 board_id_mu = stoi(it.first.substr(it.first.find("_") + 1));
844 for ( auto x : boardMapsMu["MuFilter"][it.first] )
845 {
846 for ( auto slot : slots )
847 {
848 if ( Path.find("testbeam_23") == string::npos ) s = 0;
849 else s = 3;
850 tmp = x.second.substr(0, x.second.find("_"));
851 if ( tmp =="US" ) s = 1;
852 if ( tmp=="DS" )
853 {
854 s = 2;
855 if ( Path.find("testbeam_24") != string::npos && x.second.substr(x.second.find("_")+1,1)==2) s = 3;
856 }
857 if ( slots[slot.first] == x.first )
858 {
859 MufiSystem[board_id_mu][slot.first] = s;
860 boardMaps["MuFilter"][it.first][slot.second] = slot.first;
861 }
862 }
863 }
864 }
865
866 //string o;
867 for (int i = 1 ; i < 6; i++ )
868 {
869 if (i<3)
870 {
871 offMap[Form("Veto_%iLeft",i)] = {10000 + (i-1)*1000+ 0, 8, 2};//first channel, nSiPMs, nSides
872 offMap[Form("Veto_%iRight",i)] = {10000 + (i-1)*1000+ 0, 8, 2};
873 }
874 if (i==3) offMap[Form("Veto_%iVert",i)] = {10000 + (i-1)*1000+ 6, -8, 1};// 3rd Veto plane
875 if (i<4)
876 {
877 // DS
878 offMap[Form("DS_%iLeft",i)] = {30000 + (i-1)*1000+ 59, -1, 2};// direction not known
879 offMap[Form("DS_%iRight",i)] = {30000 + (i-1)*1000+ 59, -1, 2};// direction not known
880 }
881 if (i<5)
882 {
883 offMap[Form("DS_%iVert",i)] = {30000 + (i-1)*1000+ 119, -1, 1};// direction not known
884 }
885 offMap[Form("US_%iLeft",i)] = {20000 + (i-1)*1000+ 9, -8, 2}; // from top to bottom
886 offMap[Form("US_%iRight",i)] = {20000 + (i-1)*1000+ 9, -8, 2};
887 }
888}
889/* Define calibration functions **/
890double ConvRawData::qdc_calibration( int board_id, int tofpet_id, int channel,
891 int tac, uint16_t v_coarse, uint16_t v_fine, uint16_t tf )
892 {
893 double GQDC = 1.0; // or 3.6
894 map<string, double> par= X_qdc[{board_id,tofpet_id,channel,tac}];
895 double x = v_coarse - tf;
896 double fqdc = -par["c"]*log(1+exp( par["a"]*pow((x-par["e"]),2)
897 -par["b"]*(x-par["e"]) )) + par["d"];
898 double value = (v_fine-fqdc)/GQDC;
899 return value;
900}
901double ConvRawData::qdc_chi2( int board_id, int tofpet_id, int channel, int tac, int TDC=0 )
902{
903 map<string, double> par = X_qdc[{board_id,tofpet_id,channel,tac}];
904 map<string, double> parT = X_tdc[{board_id,tofpet_id,channel,tac, TDC}];
905 return max(par["chi2Ndof"],parT["chi2Ndof"]);
906}
907double ConvRawData::qdc_sat( int board_id, int tofpet_id, int channel, int tac, uint16_t v_fine )
908{
909 map<string, double> par = X_qdc[{board_id,tofpet_id,channel,tac}];
910 return v_fine/par["d"];
911}
912double ConvRawData::time_calibration( int board_id, int tofpet_id, int channel,
913 int tac, int64_t t_coarse, uint16_t t_fine, int TDC=0 )
914{
915 map<string, double> parT = X_tdc[{board_id,tofpet_id,channel,tac, TDC}];
916 double x = t_fine;
917 double ftdc = (-parT["b"]-sqrt(pow(parT["b"],2)
918 -4*parT["a"]*(parT["c"]-x)))/(2*parT["a"]);
919 double timestamp = t_coarse+ftdc;
920 return timestamp;
921}
922tuple<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
923{
924 map<string, double> par= X_qdc[{board_id,tofpet_id,channel,tac}];
925 map<string, double> parT = X_tdc[{board_id,tofpet_id,channel,tac,TDC}];
926 double x = t_fine;
927 double ftdc = (-parT["b"]-sqrt(pow(parT["b"],2)
928 -4*parT["a"]*(parT["c"]-x)))/(2*parT["a"]); // Ettore 28/01/2022 +par['d']
929 double timestamp = t_coarse + ftdc;
930 double tf = timestamp - t_coarse;
931 x = v_coarse - tf;
932 double fqdc = - par["c"]*log(1+exp( par["a"]*pow((x-par["e"]),2)-par["b"]*(x-par["e"]) ))
933 + par["d"];
934 double value = (v_fine-fqdc)/GQDC;
935 return make_tuple(timestamp,value,max(par["chi2Ndof"],parT["chi2Ndof"]),v_fine/par["d"]);
936}
937map<double, pair<double, double> > ConvRawData::calibrationReport()
938{
939 TH1F* h = new TH1F("chi2","chi2", 1000, 0., 10000);
940 map<double, pair<double, double> > report{};
941 int TDC = 0;
942 double chi2{}, chi2T{}, key{};
943 // loop over entries of X_qdc and get map keys
944 for (auto it : X_qdc)
945 {
946 map<string, double> par= X_qdc[{it.first[0], it.first[1], it.first[2], it.first[3]}];
947 if (par["chi2Ndof"]) chi2 = par["chi2Ndof"];
948 else chi2 = -1;
949 map<string, double> parT = X_tdc[{it.first[0], it.first[1], it.first[2], it.first[3], TDC}];
950 if (parT["chi2Ndof"]) chi2T = parT["chi2Ndof"];
951 else chi2T = -1;
952 key = it.first[3] +10*it.first[2] + it.first[1]*10*100 + it.first[0]*10*100*100;
953 if (report.count(key)==0) report[key] = make_pair(chi2,chi2T);
954 }
955 for (auto it : report)
956 {
957 h->Fill(report[it.first].first);
958 h->Fill(report[it.first].second);
959 }
960 h->Draw();
961 return report;
962}
964int ConvRawData::channel_func( int tofpet_id, int tofpet_channel, int position)
965{
966 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.
967}
969void ConvRawData::read_csv(string Path)
970{
971 ifstream infile; // ifstream is the stream class to only read! from files.
972 stringstream X;
973 // reading file with xrootd protocol
974 File file;
975 XRootDStatus status;
976 StatInfo *info;
977 uint64_t offset = 0;
978 uint32_t size;
979 uint32_t bytesRead = 0;
980 string line, element;
981 double chi2_Ndof{};
982 // counter of number of elements
983 vector<int> key_vector{};
984
985 // Always read the SiPM mapping
986 int SiPM{};
987 vector<int> data_vector{};
988 map<string, map<int, vector<int>> > SiPMmap{};
989 vector<int> row{};
990 map<string, int> key { {"BM",3}, {"DS",2}, {"US",1}, {"Veto",0} };
991 struct stat buffer;
992 TString sndRoot = gSystem->Getenv("SNDSW_ROOT");
993 string sndswPath = sndRoot.Data();
994 string path_SiPMmap = Form("%s/geometry", sndswPath.c_str());
995 if (stat(path_SiPMmap.c_str(), &buffer) != 0)
996 {
997 LOG (error) << "Path "<< path_SiPMmap.c_str() << " does not exist!";
998 exit (0);
999 }
1000 for (auto sys : key)
1001 {
1002 infile.open(Form("%s/%s_SiPM_mapping.csv", path_SiPMmap.c_str(), sys.first.c_str()));
1003 X << infile.rdbuf();
1004 infile.close();
1005
1006 // Skip 1st line
1007 getline(X,line);
1008 LOG (info) << "In " << sys.first << " SiPM map file: " << line;
1009 // Read all other lines
1010 while (getline(X,line))
1011 {
1012 data_vector.clear();
1013 stringstream items(line);
1014 // first element in line is SiPM
1015 getline(items, element, ',');
1016 SiPM = stoi(element);
1017 // only taking first few elements of line
1018 while (data_vector.size()<4 && getline(items, element, ','))
1019 {
1020 data_vector.push_back(stoi(element));
1021 }
1022 SiPMmap[sys.first][SiPM] = data_vector;
1023 if (X.peek() == EOF) break;
1024 }
1025 X.str(string()); X.clear(); line.clear();
1026 size = 0; offset = 0; bytesRead = 0;
1027 for (auto channel : SiPMmap[sys.first])
1028 {
1029 row = channel.second;
1030 TofpetMap[sys.second][row.at(2)*1000+row.at(3)] = channel.first;
1031 }
1032 } // end filling SiPMmap and TofpetMap
1033
1034 // Get QDC/TDC calibration data if flag is set
1035 if (!makeCalibration) return;
1036
1037 // Get QDC calibration data
1038 if (local)
1039 {
1040 infile.open(Form("%s/qdc_cal.csv", Path.c_str()));
1041 X << infile.rdbuf();
1042 infile.close();
1043 }
1044 else
1045 {
1046 status = file.Open(Form("%s/qdc_cal.csv", Path.c_str()), OpenFlags::Read);
1047 if( !status.IsOK() )
1048 {
1049 LOG (error) << "Error opening file";
1050 exit(0);
1051 }
1052 file.Stat(false, info);
1053 size = info->GetSize();
1054 char *buff = new char[size];
1055 status = file.Read(offset, size, buff, bytesRead);
1056 X << buff;
1057 status = file.Close();
1058 delete info;
1059 delete [] buff;
1060 }
1061 LOG (info) << "Read_csv "<<Path.c_str();
1062 // Skip 1st line
1063 getline(X, line);
1064 LOG (info) << "In QDC cal file: " << line;
1065 vector<double> qdcData{};
1066 // Read all other lines
1067 while (getline(X, line))
1068 {
1069 qdcData.clear();
1070 stringstream items(line);
1071 while (getline(items, element, ','))
1072 {
1073 if(iscntrl(element[0])) break;
1074 qdcData.push_back(stof(element));
1075 }
1076 if (qdcData.size()<10) continue;
1077 key_vector.clear();
1078 key_vector = {int(qdcData[0]), int(qdcData[1]), int(qdcData[2]), int(qdcData[3])};
1079 chi2_Ndof = (qdcData[9] < 2) ? 999999. : qdcData[7]/qdcData[9];
1080 X_qdc[key_vector] = { {"a",qdcData[4]}, {"b",qdcData[5]}, {"c",qdcData[6]},
1081 {"d",qdcData[8]}, {"e",qdcData[10]}, {"chi2Ndof",chi2_Ndof} };
1082 if (X.peek() == EOF) break;
1083 }
1084 X.str(string()); X.clear(); line.clear();
1085 size = 0; offset = 0; bytesRead = 0;
1086 // Get TDC calibration data
1087 if (local)
1088 {
1089 infile.open(Form("%s/tdc_cal.csv", Path.c_str()));
1090 X << infile.rdbuf();
1091 infile.close();
1092 }
1093 else
1094 {
1095 status = file.Open(Form("%s/tdc_cal.csv", Path.c_str()), OpenFlags::Read);
1096 if( !status.IsOK() )
1097 {
1098 LOG (error) << "Error opening file";
1099 exit(0);
1100 }
1101 file.Stat(false, info);
1102 size = info->GetSize();
1103 char *buff = new char[size];
1104 status = file.Read(offset, size, buff, bytesRead);
1105 X << buff;
1106 status = file.Close();
1107 delete info;
1108 delete [] buff;
1109 }
1110 // Skip 1st line
1111 getline(X, line);
1112 LOG (info) << "In TDC cal file: " << line;
1113 vector<double> tdcData{};
1114 // Read all other lines
1115 while (getline(X, line))
1116 {
1117 tdcData.clear();
1118 stringstream items(line);
1119 if (line.length()<5){continue;}
1120 while (getline(items, element, ','))
1121 {
1122 if(iscntrl(element[0])) break;
1123 tdcData.push_back(stof(element));
1124 }
1125 if (tdcData.size()<9) continue;
1126 key_vector.clear();
1127 key_vector = {int(tdcData[0]), int(tdcData[1]), int(tdcData[2]), int(tdcData[3]), int(tdcData[4])};
1128 chi2_Ndof = (tdcData[10] < 2) ? 999999. : tdcData[8]/tdcData[10];
1129 X_tdc[key_vector] = { {"a",tdcData[5]}, {"b",tdcData[6]}, {"c",tdcData[7]},
1130 {"d",tdcData[9]}, {"chi2Ndof",chi2_Ndof} };
1131 if (X.peek() == EOF) break;
1132 }
1133 X.str(string()); X.clear(); line.clear();
1134 size = 0; offset = 0; bytesRead = 0;
1135}
1136void ConvRawData::debugMapping(string brd, int tofpetID, int tofpetChannel)
1137{
1138 int Key{}, SiPMChannel{}, n_SiPMs{}, n_Sides{}, Direction{}, DetID{}, SiPM_number{};
1139 string Tmp;
1140 int brdID{}, System{};
1141 Key = (tofpetID%2)*1000 + tofpetChannel;
1142 Tmp = boardMapsMu["MuFilter"][brd][slots[tofpetID]];
1143 brdID = stoi(brd.substr(brd.find("_")+1));
1144 System = MufiSystem[brdID][tofpetID];
1145 SiPMChannel = TofpetMap[System][Key]-1;
1146 n_SiPMs = abs(offMap[Tmp][1]);
1147 n_Sides = abs(offMap[Tmp][2]);
1148 Direction = int(offMap[Tmp][1]/n_SiPMs);
1149 DetID = offMap[Tmp][0] + Direction*int(SiPMChannel/n_SiPMs);
1150 SiPM_number = SiPMChannel%(n_SiPMs);
1151 cout << "SiPM channel " << SiPMChannel << " nSiPM " << n_SiPMs << " nSides " << n_Sides
1152 << " detID " << DetID << " SiPM number " << SiPM_number << endl;
1153}
1154
1155ClassImp(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:901
ConvRawData::fEventHeader
FairEventHeader * fEventHeader
Definition ConvRawData.h:104
ConvRawData::debugMapping
void debugMapping(std::string board, int tofpetID, int tofpetChannel)
Definition ConvRawData.cxx:1136
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:890
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:912
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:801
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:754
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:922
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:907
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:964
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:969
ConvRawData::makeCalibration
int makeCalibration
Definition ConvRawData.h:94
ConvRawData::calibrationReport
std::map< double, std::pair< double, double > > calibrationReport()
Definition ConvRawData.cxx:937
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