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