ATLAS Offline Software
Loading...
Searching...
No Matches
Calib.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5//***********************************************************************************************
6// PixelCalibration - Program to run the pixel charge calibration (not IBL)
7// ---------------------------------------
8// begin : 01 11 2023
9// email : sergi.rodriguez@cern.ch
10//***********************************************************************************************
11
12
13#include "Calib.h"
14#include <format>
15
16bool Calib::totFitting (const pix::PixelMapping &pm, const std::string &inTotFile, std::map<unsigned int , std::vector<std::unique_ptr<CalibFrontEndInfo>> > &map_info){
17
18 if (inTotFile.empty()) return false;
19
20 TFile totFile(inTotFile.c_str(),"READ");
21 if (not totFile.IsOpen()) {
22 printf("Error - File %s could not be opened.\n",inTotFile.c_str());
23 totFile.Close();
24 return false;
25 } else {
26 printf("File %s opened.\n",inTotFile.c_str());
27 printf("Running TOT calibration...\n");
28 }
29
30 // Creating here the histograms with this scope ownership
31 std::vector< std::vector< std::unique_ptr<TH1F> > > histogramsTOT;
32 std::vector< std::vector< std::unique_ptr<TH1F> > > histogramsTOTsig;
33 for(unsigned int FE = 0; FE < m_nFE; FE++){
34
35 histogramsTOT.push_back( std::vector< std::unique_ptr<TH1F> >() );
36 histogramsTOTsig.push_back( std::vector< std::unique_ptr<TH1F> >() );
37
38 //Here we combine long and ganged pixels
39 for(unsigned int pixel=0; pixel<2; pixel++){
40 std::string title = "FE"+std::to_string(FE)+"_pixType"+std::to_string(pixel);
41 std::unique_ptr<TH1F> htot = std::make_unique<TH1F>((title+"_tot").c_str(), (title+"_tot").c_str(), m_totnbins, m_totLo, m_totHi);
42 htot->SetDirectory(0);
43 std::unique_ptr<TH1F> htotsig = std::make_unique<TH1F>((title+"_totsig").c_str(), (title+"_totsig").c_str(), m_totsigNBins, m_totsigLo, m_totsigHi);
44 htotsig->SetDirectory(0);
45 //cppcheck-suppress containerOutOfBounds
46 histogramsTOT.at(FE).push_back(std::move(htot));
47 //cppcheck-suppress containerOutOfBounds
48 histogramsTOTsig.at(FE).push_back(std::move(htotsig));
49 }
50 }
51
52 // Start looping over the ROD, MOD and charges
53 TIter rodItr=getRodIterator(totFile);
54 TKey* rodKey;
55
56 while ((rodKey=static_cast<TKey*>(rodItr()))) {
57 TDirectoryFile* rodDir = static_cast<TDirectoryFile*>(rodKey->ReadObj());
58 TIter modItr = getModuleIterator(rodDir);
59 const std::string rodName(rodKey->GetName());
60 printf("%s\n",rodName.c_str());
61 TKey* modKey;
62 while ((modKey=static_cast<TKey*>(modItr()))) {
63 const std::string modName(modKey->GetName());
64
65 if (not moduleInPart(modName)) continue;
66 if (not pm.contains(modName)) continue;
67
68 if( m_runOneMOD and (modName == m_testMOD)){
69 continue;
70 }
71
72 printf(" -> %s\n",modName.c_str());
73
74 //creates arrays for the Tgraph
75 std::array<std::array<float, m_ncharge>, m_nFE> totArrI{};
76 std::array<std::array<float, m_ncharge>, m_nFE> totErrArrI{};
77 std::array<std::array<float, m_ncharge>, m_nFE> totSigArrI{};
78 std::array<std::array<float, m_ncharge>, m_nFE> totSigErrArrI{};
79 std::array<std::array<float, m_ncharge>, m_nFE> totLongArrI{};
80 std::array<std::array<float, m_ncharge>, m_nFE> totErrLongArrI{};
81
82 // loop over charges
83 static const std::string totMeanStr{"TOT_MEAN"};
84 static const std::string totSigmaStr{"TOT_SIGMA"};
85 for (int c=0; c<m_ncharge; ++c) {
86
87 // Get TH2 for a given charge
88 std::unique_ptr<TH2F> h2dTOTmean(get2DHistogramFromPath(rodDir,modName, totMeanStr, c));
89 std::unique_ptr<TH2F> h2dTOTsig(get2DHistogramFromPath(rodDir,modName, totSigmaStr, c));
90 if(!h2dTOTmean or !h2dTOTsig) {
91 return false;
92 }
93 h2dTOTmean->SetDirectory(0);
94 h2dTOTsig->SetDirectory(0);
95 // loop over pixels
96 for (unsigned int ieta = 0; ieta < m_etaBins; ieta++) {
97 for (unsigned int iphi = 0; iphi < m_phiBins; iphi++) {
98 float totmean = h2dTOTmean->GetBinContent(ieta + 1, iphi + 1);
99 float totsig = h2dTOTsig ->GetBinContent(ieta + 1, iphi + 1);
100
101 if (totmean<0.1) {
102 continue;
103 }
104
105 int FE = chipId(iphi, ieta);
106 int pixel= pixelType(iphi, ieta, true);
107
108 if(FE<0){
109 return false;
110 }
111
112 histogramsTOT.at(FE).at(pixel)->Fill(totmean);
113 histogramsTOTsig.at(FE).at(pixel)->Fill(totsig);
114
115 }
116 }
117
118 // free memory
119 h2dTOTmean.reset();
120 h2dTOTsig.reset();
121
122 //filling arrays
123 for(unsigned int FE = 0; FE < m_nFE; FE++){
124 for(unsigned int pixel = 0; pixel <2; pixel++){
125
126 if(pixel == 0){
127 totArrI.at(FE).at(c) = histogramsTOT.at(FE).at(pixel)->GetMean();
128 totErrArrI.at(FE).at(c) = histogramsTOT.at(FE).at(pixel)->GetMeanError();
129 totSigArrI.at(FE).at(c) = std::sqrt(std::pow(histogramsTOTsig.at(FE).at(pixel)->GetMean() ,2)+std::pow(histogramsTOT.at(FE).at(pixel)->GetRMS() ,2));
130 totSigErrArrI.at(FE).at(c) = std::sqrt(std::pow(histogramsTOTsig.at(FE).at(pixel)->GetMeanError(),2)+std::pow(histogramsTOT.at(FE).at(pixel)->GetRMSError(),2));
131
132 if(totSigErrArrI.at(FE).at(c) > 1.0){
133 totArrI.at(FE).at(c) = 0.0;
134 }
135
136 }
137 else{
138 totLongArrI.at(FE).at(c) = histogramsTOT.at(FE).at(pixel)->GetMean();
139 totErrLongArrI.at(FE).at(c) = histogramsTOT.at(FE).at(pixel)->GetMeanError();
140
141 if(totErrLongArrI.at(FE).at(c) > 1.0){
142 totLongArrI.at(FE).at(c) = 0.0;
143 }
144
145 }
146
147 //reset histogram for next iteration
148 histogramsTOT.at(FE).at(pixel)->Reset("ICESM");
149 histogramsTOTsig.at(FE).at(pixel)->Reset("ICESM");
150 }
151 }
152 } // End of charge.
153
154 // loop over FE and create a graph for fitting
155 for(unsigned int FE = 0; FE < m_nFE; FE++) {
156
157 std::string subdir(((FE < 10) ? "FE0" : "FE") +std::to_string(FE));
158
159 std::vector<float> v_Q;
160 std::vector<float> v_Qerr;
161 std::vector<float> v_TOT;
162 std::vector<float> v_TOTerr;
163 std::vector<float> v_TOTsig;
164 std::vector<float> v_TOTsigerr;
165 std::vector<float> v_TOTlong;
166 std::vector<float> v_TOTlongerr;
167
168 std::copy(std::begin(m_chargeArr) , std::end(m_chargeArr) , std::back_inserter(v_Q) );
169 std::copy(std::begin(m_chargeErrArr) , std::end(m_chargeErrArr) , std::back_inserter(v_Qerr) );
170 std::copy(std::begin(totArrI[FE]) , std::end(totArrI[FE]) , std::back_inserter(v_TOT) );
171 std::copy(std::begin(totErrArrI[FE]) , std::end(totErrArrI[FE]) , std::back_inserter(v_TOTerr) );
172 std::copy(std::begin(totSigArrI[FE]) , std::end(totSigArrI[FE]) , std::back_inserter(v_TOTsig) );
173 std::copy(std::begin(totSigErrArrI[FE]) , std::end(totSigErrArrI[FE]) , std::back_inserter(v_TOTsigerr) );
174 std::copy(std::begin(totLongArrI[FE]) , std::end(totLongArrI[FE]) , std::back_inserter(v_TOTlong) );
175 std::copy(std::begin(totErrLongArrI[FE]) , std::end(totErrLongArrI[FE]) , std::back_inserter(v_TOTlongerr) );
176
177 std::vector<float> pixNormalParams;
178 std::vector<float> pixNormalParamsQuality;
179 std::vector<float> pixSigParams;
180 std::vector<float> pixSigParamsQuality;
181 std::vector<float> pixLongParams;
182 std::vector<float> pixLongParamsQuality;
183
184 //For normal pixels and sig
185 uint8_t n_fit = 0;
186 do{
187 int vecsize = v_Q.size();
188
189 std::unique_ptr<TGraphErrors> graphnormal = std::make_unique<TGraphErrors>(vecsize, &v_Q.at(0), &v_TOT.at(0) , &v_Qerr.at(0), &v_TOTerr.at(0) );
190 std::unique_ptr<TGraphErrors> graphsig = std::make_unique<TGraphErrors>(vecsize, &v_Q.at(0), &v_TOTsig.at(0) , &v_Qerr.at(0), &v_TOTsigerr.at(0) );
191
192 std::unique_ptr<TF1> functnormal = std::make_unique<TF1>("normal" ,new funcTot , m_chargeArr[m_qthresh]-100, m_chargeArr[m_ncharge-1]+100, 3);
193 std::unique_ptr<TF1> functnormalsig = std::make_unique<TF1>("normal_sig",new funcDisp, m_chargeArr[m_qthresh]-100, m_chargeArr[m_ncharge-1]+100, 2);
194
195 graphnormal->Fit(functnormal.get() ,"MRQ");
196 graphsig ->Fit(functnormalsig.get(),"MRQ");
197
198 pixNormalParams = getParams(functnormal.get() ,3 );
199 pixSigParams = getParams(functnormalsig.get(),2 );
200
201 pixNormalParamsQuality = getParams_quality(functnormal.get() );
202 pixSigParamsQuality = getParams_quality(functnormalsig.get());
203
204 if(m_savefile){
205
206 m_wFile->cd();
207 if( !m_wFile->Get((rodName+"/"+modName+"/TOTfits/"+subdir).c_str()) ){
208 m_wFile->mkdir((rodName+"/"+modName+"/TOTfits/"+subdir).c_str(),rodName.c_str());
209 }
210
211 m_wFile->cd((rodName+"/"+modName+"/TOTfits/"+subdir).c_str());
212
213 graphTitles(graphnormal, std::string(modName+" - "+subdir+" - normal pixels: Fit: "+std::to_string(n_fit)).c_str(), "TOT");
214 graphTitles(graphsig , std::string(modName+" - "+subdir+" - normal pixels: Fit: "+std::to_string(n_fit)).c_str(), "Charge smearing");
215
216 graphnormal->Write(std::string("normal_fit_"+std::to_string(n_fit)).c_str(), TObject::kWriteDelete);
217 graphsig->Write(std::string("smearing_fit_"+std::to_string(n_fit)).c_str(), TObject::kWriteDelete);
218 n_fit++;
219 }
220
221 functnormal.reset();
222 functnormalsig.reset();
223
224 graphnormal.reset();
225 graphsig.reset();
226
227 }while(reFit_normalPix(pixNormalParams, v_Q, v_Qerr, v_TOT, v_TOTerr, v_TOTsig, v_TOTsigerr, FE ) );
228
229
230 // Since we have modified the vector size we need to clear it and refill it for the long and gange pixels
231 v_Q.clear();
232 v_Qerr.clear();
233
234 std::copy(std::begin(m_chargeArr) , std::end(m_chargeArr) , std::back_inserter(v_Q) );
235 std::copy(std::begin(m_chargeErrArr) , std::end(m_chargeErrArr) , std::back_inserter(v_Qerr) );
236
237
238 //For long and ganged pixels (combined)
239 n_fit = 0;
240 do{
241
242 //size will be modified in the condition..
243 int vecsize = v_Q.size();
244
245 std::unique_ptr<TGraphErrors> graflong = std::make_unique<TGraphErrors>(vecsize, &v_Q.at(0), &v_TOTlong.at(0), &v_Qerr.at(0), &v_TOTlongerr.at(0) );
246
247 std::unique_ptr<TF1> functlong = std::make_unique<TF1>("long" ,new funcTot , m_chargeArr[m_qthresh]-100, m_chargeArr[m_ncharge-1]+100, 3);
248
249 graflong ->Fit(functlong.get() ,"MRQ");
250
251 pixLongParams = getParams(functlong.get() ,3 );
252 pixLongParamsQuality = getParams_quality(functlong.get() );
253
254 if(m_savefile){
255
256 m_wFile->cd();
257 if( !m_wFile->Get((rodName+"/"+modName+"/TOTfits/"+subdir).c_str()) ){
258 m_wFile->mkdir((rodName+"/"+modName+"/TOTfits/"+subdir).c_str(),rodName.c_str());
259 }
260
261 m_wFile->cd((rodName+"/"+modName+"/TOTfits/"+subdir).c_str());
262
263 graphTitles(graflong, std::string(modName+" - "+subdir+" - long+ganged pixels: Fit: "+std::to_string(n_fit)).c_str(), "TOT");
264
265 graflong->Write(std::string("long_ganged_fit_"+std::to_string(n_fit)).c_str(), TObject::kWriteDelete);
266 n_fit++;
267 }
268
269 //delete the TF1
270 functlong.reset();
271
272 //delete the graphs
273 graflong.reset();
274
275 // no need to loop here.. leaving to improve it in the future - might need refitting
276 }while( false);
277
278
279 // Find the module it belong to
280 int modID = pm.getID(modName);
281 auto itr = map_info.find( modID );
282
283 if (itr != map_info.end()) {
284 (itr->second).at(FE)->set_NormalParams( pixNormalParams);
285 (itr->second).at(FE)->set_LongParams ( pixLongParams );
286 (itr->second).at(FE)->set_SigParams ( pixSigParams );
287
288 (itr->second).at(FE)->set_times_fitted ( n_fit );
289
290 (itr->second).at(FE)->set_NormalParamsQuality( pixNormalParamsQuality);
291 (itr->second).at(FE)->set_LongParamsQuality ( pixLongParamsQuality );
292 (itr->second).at(FE)->set_SigParamsQuality ( pixSigParamsQuality );
293 }
294 else{
295 printf("Error - Module not found in fitting step... Skipping.\n");
296 return false;
297 }
298
299 } // End of FE loop
300 } // End of MOD
301 } // End of ROD
302
303 // remove from memory
304 for(unsigned int FE = 0; FE < m_nFE; FE++) {
305 for(unsigned int pixel=0; pixel<2; pixel++) {
306 histogramsTOT.at(FE).at(pixel).reset();
307 histogramsTOTsig.at(FE).at(pixel).reset();
308 }
309 }
310
311 totFile.Close();
312 return true;
313}
314
315
316
317bool Calib::reFit_normalPix(std::vector<float> &params, std::vector<float> &q, std::vector<float> &qerr, std::vector<float> &tot, std::vector<float> &toterr, std::vector<float> &sig, std::vector<float> &sigerr, const unsigned int fe){
318
319 float vecFit_size = q.size()+1;
320 // float vecFit_size = q.size() - m_qthresh;
321 float stopFit = m_ncharge/2.0;
322 // float stopFit = (m_ncharge - m_qthresh)/2.0;
323 if(vecFit_size < stopFit) {
324
325 // Default values for the fit
326 params.at(0) = 0;
327 params.at(1) = -28284.3;
328 params.at(2) = 0;
329
330 printf("reFit_normalPix: Refitting skipped. Not enough points to fit.\n");
331
332 return false;
333 }
334
335 float parAI0 = params.at(0);
336 float parEI0 = params.at(1);
337 float parCI0 = params.at(2);
338
339 std::vector<float> v_discrepancy;
340
341 for(unsigned int i = 0; i < q.size(); i++){
342 float discrepancy = std::abs( 1 - ( (parAI0 * parEI0 - parCI0 * tot.at(i)) / (tot.at(i) - parAI0) ) / q.at(i) );
343
344 if( i < m_qthresh ){
345 discrepancy = 0.0;
346 }
347 v_discrepancy.push_back(discrepancy);
348 }
349
350 auto itr_max = std::max_element(v_discrepancy.begin(),v_discrepancy.end());
351
352 if(*itr_max > m_chi_error){
353
354 size_t n_max = std::distance(v_discrepancy.begin(), itr_max);
355 printf("FE %02u Refitted, removing charge %5.0f with chi_error %7.5f\n", fe ,q.at(n_max),*itr_max);
356 q.erase(q.begin()+n_max);
357 qerr.erase(qerr.begin()+n_max);
358 tot.erase(tot.begin()+n_max);
359 toterr.erase(toterr.begin()+n_max);
360 sig.erase(sig.begin()+n_max);
361 sigerr.erase(sigerr.begin()+n_max);
362
363 return true;
364 }
365
366 return false;
367
368}
369
370
371bool Calib::fillTiming(const pix::PixelMapping &pm, const std::string &inTimFile, std::map<unsigned int, std::vector<std::unique_ptr<CalibFrontEndInfo>> > &map_info) {
372
373 if (inTimFile.empty()) return false;
374
375 TFile timFile(inTimFile.c_str(),"READ");
376 if (not timFile.IsOpen()) {
377 printf("Error - File %s could not be opened.\n",inTimFile.c_str());
378 return false;
379 } else {
380 printf("File %s opened.\n",inTimFile.c_str());
381 printf("Running timming calibration...\n");
382 }
383
384 // Creating here the histograms with this scope owner
385 std::vector< std::vector< std::unique_ptr<TH1F> > > histogramsTIM;
386 for(unsigned int FE = 0; FE < m_nFE; FE++){
387
388 histogramsTIM.push_back( std::vector< std::unique_ptr<TH1F> >() );
389
390 for(unsigned int pixel=0; pixel<3; pixel++){
391 std::string title = "FE"+std::to_string(FE)+"_pixType"+std::to_string(pixel);
392 std::unique_ptr<TH1F> h = std::make_unique<TH1F>((title+"_thr").c_str(), (title+"_thr").c_str(), m_timnbins, m_timLo, m_timHi);
393 h->SetDirectory(0);
394 //cppcheck-suppress containerOutOfBounds
395 histogramsTIM.at(FE).push_back(std::move(h));
396 }
397 }
398
399 //Start looping over the RODs
400 TIter rodItr = getRodIterator(timFile);
401 TKey* rodKey;
402 while ((rodKey=static_cast<TKey*>(rodItr()))) {
403 const std::string rodName(rodKey->GetName());
404 TDirectoryFile* rodDir = static_cast<TDirectoryFile*>(rodKey->ReadObj());
405 TKey* modKey;
406 TIter modItr=getModuleIterator(rodDir);
407
408 // Looping over the MODs of each ROD
409 static const std::string sCurveMeanStr{"SCURVE_MEAN"};
410 while ((modKey=static_cast<TKey*>(modItr()))) {
411 std::string modName(modKey->GetName());
412
413 if ( not moduleInPart(modName)){
414 continue;
415 }
416 if ( not pm.contains(modName)){
417 printf("Error - Module %s not found in the PixelMapping tool\n",modName.c_str());
418 continue;
419 }
420
421 if( m_runOneMOD and (modName == m_testMOD)){
422 continue;
423 }
424
425 std::unique_ptr<TH2F> h2dTim(get2DHistogramFromPath(rodDir,modName, sCurveMeanStr));
426 h2dTim->SetDirectory(0);
427
428 for (unsigned int ieta = 0; ieta < m_etaBins; ieta++) {
429 for (unsigned int iphi = 0; iphi < m_phiBins; iphi++) {
430
431 float tim = h2dTim->GetBinContent(ieta + 1, iphi + 1);
432
433 if (tim<0.5) {
434 continue;
435 }
436
437 int FE = chipId(iphi, ieta);
438 int pixel= pixelType(iphi, ieta);
439
440 if(FE<0){
441 return false;
442 }
443
444 histogramsTIM.at(FE).at(pixel)->Fill(tim);
445 }
446 }
447
448 // Freeing memory of TH2F
449 h2dTim.reset();
450
451 int modID = pm.getID(modName);
452 auto itr = map_info.find( modID );
453 if (itr == map_info.end()) {
454 printf("Calib::fillTiming: ERROR - Mod ID= %16s not found. Creating it -----> Inform Pixel Offline Software Experts... \n",modName.c_str());
455
456 map_info[modID] = std::vector<std::unique_ptr<CalibFrontEndInfo>> ();
457
458 for(unsigned int FE = 0; FE < m_nFE; FE++){
459
460 map_info[modID].push_back( std::unique_ptr<CalibFrontEndInfo>() );
461 std::unique_ptr<CalibFrontEndInfo> p = std::make_unique<CalibFrontEndInfo>(modID,FE,modName,std::string(rodKey->GetName()));
462 map_info[modID].at(FE) = std::move(p);
463
464 for(unsigned int pixel=0; pixel<3; pixel++){
465
466 // Saving information for the calibration
467 int tim_mean = histogramsTIM.at(FE).at(pixel)->GetMean();
468
469 // Reset histograms for next front end
470 histogramsTIM.at(FE).at(pixel)->Reset("ICESM");
471
472 if(pixel == 0){ // normal
473 map_info[modID].at(FE)->set_NormalIntime(tim_mean);
474 }
475 else if(pixel == 1){ // long
476 map_info[modID].at(FE)->set_LongIntime(tim_mean);
477 }
478 else if(pixel == 2){ // ganged
479 map_info[modID].at(FE)->set_GangedIntime(tim_mean);
480 }
481 else{
482 printf("Error - Bad pixel in Calib::fillTiming\n");
483 return false;
484 }
485 }
486 }
487
488 }
489 else{
490 for(unsigned int FE = 0; FE < m_nFE; FE++){
491 std::string subdir(((FE < 10) ? "FE0" : "FE") +std::to_string(FE));
492 for(unsigned int pixel=0; pixel<3; pixel++){
493
494 // Saving information for the calibration
495 int tim_mean = histogramsTIM.at(FE).at(pixel)->GetMean();
496
497 if(m_savefile){
498
499 m_wFile->cd();
500 if( !m_wFile->Get((rodName+"/"+modName+"/Thresholds/"+subdir).c_str()) ){
501 m_wFile->mkdir((rodName+"/"+modName+"/Thresholds/"+subdir).c_str(),rodName.c_str());
502 }
503
504 m_wFile->cd((rodName+"/"+modName+"/Thresholds/"+subdir).c_str());
505
506 histogramsTIM.at(FE).at(pixel)->SetTitle("Intime;Pixel intime;Counts");
507
508 std::string type = "";
509 if(pixel == 0) type = "normal";
510 else if(pixel == 1) type = "long";
511 else if(pixel == 2) type = "ganged";
512
513 histogramsTIM.at(FE).at(pixel)->Write(std::string("intime_"+type).c_str(), TObject::kWriteDelete);
514 }
515
516 // Reset histograms for next front end
517 histogramsTIM.at(FE).at(pixel)->Reset("ICESM");
518
519 if(pixel == 0){ // normal
520 (itr->second).at(FE)->set_NormalIntime(tim_mean);
521 }
522 else if(pixel == 1){ // long
523 (itr->second).at(FE)->set_LongIntime(tim_mean);
524 }
525 else if(pixel == 2){ // ganged
526 (itr->second).at(FE)->set_GangedIntime(tim_mean);
527 }
528 else{
529 printf("Error - Bad pixel in Calib::fillTiming\n");
530 return false;
531 }
532 } // End of pixel type loop
533 } // End of FE loop
534 }
535 } // End of MOD loop
536 }// End of ROD loop
537
538 for(unsigned int FE = 0; FE < m_nFE; FE++) {
539 for(unsigned int pixel=0; pixel<3; pixel++) {
540 histogramsTIM.at(FE).at(pixel).reset();
541 }
542 }
543
544 //Closing file - not needed anymore
545 timFile.Close();
546 printf("DONE with threshold calibration.\n");
547 return true;
548}
549
550
551
552bool Calib::fillThresholds(const pix::PixelMapping &pm, const std::string &inThrFile, std::map<unsigned int, std::vector<std::unique_ptr<CalibFrontEndInfo>> > &map_info) {
553
554 if (inThrFile.empty()) return false;
555
556 TFile riThrFile(inThrFile.c_str(),"READ");
557 if (not riThrFile.IsOpen()) {
558 printf("Error - File %s could not be opened.\n",inThrFile.c_str());
559 return false;
560 } else {
561 printf("File %s opened.\n",inThrFile.c_str());
562 printf("Running threshold calibration...\n");
563 }
564
565 // Creating here the histograms with this scope owner
566 std::vector< std::vector< std::unique_ptr<TH1F> > > histogramsTHR;
567 std::vector< std::vector< std::unique_ptr<TH1F> > > histogramsSIG;
568 for(unsigned int FE = 0; FE < m_nFE; FE++){
569
570 histogramsTHR.push_back( std::vector< std::unique_ptr<TH1F> >() );
571 histogramsSIG.push_back( std::vector< std::unique_ptr<TH1F> >() );
572
573 for(unsigned int pixel=0; pixel<3; pixel++){
574 std::string title = "FE"+std::to_string(FE)+"_pixType"+std::to_string(pixel);
575 std::unique_ptr<TH1F> hthr = std::make_unique<TH1F>((title+"_thr").c_str(), (title+"_thr").c_str(), m_thrnbins, m_thrLo, m_thrHi);
576 hthr->SetDirectory(0);
577 std::unique_ptr<TH1F> hsig = std::make_unique<TH1F>((title+"_sig").c_str(), (title+"_sig").c_str(), m_thrnbins, m_sigLo, m_sigHi);
578 hsig->SetDirectory(0);
579 //cppcheck-suppress containerOutOfBounds
580 histogramsTHR.at(FE).push_back(std::move(hthr));
581 //cppcheck-suppress containerOutOfBounds
582 histogramsSIG.at(FE).push_back(std::move(hsig));
583 }
584 }
585
586 //Will start looping over the RODs
587 TIter rodItr = getRodIterator(riThrFile);
588 TKey* rodKey;
589 static const std::string sCurveMeanStr{"SCURVE_MEAN"};
590 static const std::string sCurveSigmaStr{"SCURVE_SIGMA"};
591 while ((rodKey=static_cast<TKey*>(rodItr()))) {
592 const std::string rodName(rodKey->GetName());
593 TDirectoryFile* rodDir = (TDirectoryFile*)rodKey->ReadObj();
594 TKey* modKey;
595 TIter modItr=getModuleIterator(rodDir);
596
597 // Looping over the MODs of each ROD
598 while ((modKey=static_cast<TKey*>(modItr()))) {
599 std::string modName(modKey->GetName());
600
601 if ( not moduleInPart(modName)){
602 continue;
603 }
604 if ( not pm.contains(modName)){
605 printf("Error - Module %s not found in the PixelMapping tool\n",modName.c_str());
606 continue;
607 }
608
609 if( m_runOneMOD and (modName == m_testMOD)){
610 continue;
611 }
612
613 // pixel discriminator threshold
614 std::unique_ptr<TH2F> h2dThr(get2DHistogramFromPath(rodDir,modName, sCurveMeanStr));
615 h2dThr->SetDirectory(0);
616
617 // Getting histogram for noise
618 std::unique_ptr<TH2F>h2dSig(get2DHistogramFromPath(rodDir,modName, sCurveSigmaStr));
619 h2dSig->SetDirectory(0);
620
621 for (unsigned int ieta = 0; ieta < m_etaBins; ieta++) {
622 for (unsigned int iphi = 0; iphi < m_phiBins; iphi++) {
623
624 float thr = h2dThr->GetBinContent(ieta + 1, iphi + 1);
625 float sig = h2dSig->GetBinContent(ieta + 1, iphi + 1);
626
627 if (thr == 0 || thr > 10000 || sig == 0 || sig > 1000) {
628 continue;
629 }
630
631 int FE = chipId(iphi, ieta);
632 int pixel= pixelType(iphi, ieta);
633
634 if(FE<0){
635 return false;
636 }
637
638 histogramsTHR.at(FE).at(pixel)->Fill(thr);
639 histogramsSIG.at(FE).at(pixel)->Fill(sig);
640
641 }
642 }
643
644 // Freeing memory of TH2F
645 h2dThr.reset();
646 h2dSig.reset();
647
648 int modID = pm.getID(modName);
649 auto itr = map_info.find( modID );
650
651 // Map should be empty and therefore we need to create the key - if the key is repeated then it will throw an error
652 if (itr == map_info.end()) {
653
654 map_info[modID] = std::vector<std::unique_ptr<CalibFrontEndInfo>> ();
655
656 for(unsigned int FE = 0; FE < m_nFE; FE++){
657 std::string subdir(((FE < 10) ? "FE0" : "FE") +std::to_string(FE));
658 map_info[modID].push_back( std::unique_ptr<CalibFrontEndInfo>() );
659 std::unique_ptr<CalibFrontEndInfo> p = std::make_unique<CalibFrontEndInfo>(modID,FE,modName,std::string(rodKey->GetName()));
660 map_info[modID].at(FE) = std::move(p);
661
662 for(unsigned int pixel=0; pixel<3; pixel++){
663
664 // Saving information for the calibration
665 int thr_mean = histogramsTHR.at(FE).at(pixel)->GetMean();
666 int thr_rms = histogramsTHR.at(FE).at(pixel)->GetRMS();
667 int sig_mean = histogramsSIG.at(FE).at(pixel)->GetMean();
668
669 if(m_savefile){
670
671 m_wFile->cd();
672 if( !m_wFile->Get((rodName+"/"+modName+"/Thresholds/"+subdir).c_str()) ){
673 m_wFile->mkdir((rodName+"/"+modName+"/Thresholds/"+subdir).c_str(),rodName.c_str());
674 }
675
676 m_wFile->cd((rodName+"/"+modName+"/Thresholds/"+subdir).c_str());
677
678 histogramsTHR.at(FE).at(pixel)->SetTitle("Threshold;Pixel threshold;Counts");
679 histogramsSIG.at(FE).at(pixel)->SetTitle("Sigma;Pixel sigma;Counts");
680
681 std::string type = "";
682 if(pixel == 0) type = "normal";
683 else if(pixel == 1) type = "long";
684 else if(pixel == 2) type = "ganged";
685
686 histogramsTHR.at(FE).at(pixel)->Write(std::string("thres_"+type).c_str(), TObject::kWriteDelete);
687 histogramsSIG.at(FE).at(pixel)->Write(std::string("sigma_"+type).c_str(), TObject::kWriteDelete);
688 }
689
690 // Reset histograms for next front end
691 histogramsTHR.at(FE).at(pixel)->Reset("ICESM");
692 histogramsSIG.at(FE).at(pixel)->Reset("ICESM");
693
694 if(pixel == 0){ // normal
695 map_info[modID].at(FE)->set_NormalThreshold(thr_mean);
696 map_info[modID].at(FE)->set_NormalRms(thr_rms);
697 map_info[modID].at(FE)->set_NormalNoise(sig_mean);
698 }
699 else if(pixel == 1){ // long
700 map_info[modID].at(FE)->set_LongThreshold(thr_mean);
701 map_info[modID].at(FE)->set_LongRms(thr_rms);
702 map_info[modID].at(FE)->set_LongNoise(sig_mean);
703 }
704 else if(pixel == 2){ // ganged
705 map_info[modID].at(FE)->set_GangedThreshold(thr_mean);
706 map_info[modID].at(FE)->set_GangedRms(thr_rms);
707 map_info[modID].at(FE)->set_GangedNoise(sig_mean);
708 }
709 else{
710 printf("Calib::fillThresholds: ERROR - Bad pixel in Calib::fillThresholds\n");
711 return false;
712 }
713 }
714 }
715
716 }
717 else{
718 printf("Calib::fillThresholds: ERROR - REPEATED MOD ID: %s! Contact Offline team\n",modName.c_str());
719 return false;
720 }
721 } // End of MOD loop
722 } // End of ROD loop
723
724 for(unsigned int FE = 0; FE < m_nFE; FE++) {
725 for(unsigned int pixel=0; pixel<3; pixel++) {
726 histogramsTHR.at(FE).at(pixel).reset();
727 histogramsSIG.at(FE).at(pixel).reset();
728 }
729 }
730
731 //Closing file - not needed anymore
732 riThrFile.Close();
733 printf("DONE with threshold calibration.\n");
734 return true;
735}
736
737
738int Calib::chipId(int iphi, int ieta) {
739 int circ = -1;
740 if (iphi < 160) {
741 circ = (int)(ieta / 18);
742 } else {
743 circ = 15 - (int)(ieta / 18);
744 } // FE15, FE14, ... FE8
745
746 if (circ>15){
747 printf("Error - FE id error: %d, setting it to -1",circ);
748 circ = -1;//error
749 }
750 return circ;
751}
752
753
754int Calib::pixelType(int iphi, int ieta, bool isForTOT) {
755
756 // normal pixels ( by default )
757 int pixtype = 0;
758
759 // define long pixels
760 if (ieta % 18 == 0 || ieta % 18 == 17) {
761 pixtype = 1;
762 }
763 // define ganged pixels
764 if (iphi > 152 && iphi < 160 && iphi % 2 == 1) {
765 pixtype = 2;
766 }
767 if (iphi > 159 && iphi < 167 && iphi % 2 == 0) {
768 pixtype = 2;
769 }
770
771 if(isForTOT and pixtype == 2 ) pixtype=1;
772 return pixtype;
773}
774
775TIter Calib::getRodIterator(const TFile & inputFile) {
776 TDirectoryFile* scanDir = static_cast<TDirectoryFile*>((static_cast<TKey*>(inputFile.GetListOfKeys()->First()))->ReadObj());
777 TList* rodKeyList = static_cast<TList*>(scanDir->GetListOfKeys());
778 return TIter(rodKeyList);
779}
780
781TIter Calib::getModuleIterator( TDirectoryFile* rodDir) {
782 TList* modKeyList = static_cast<TList*>(rodDir->GetListOfKeys());
783 return TIter(modKeyList);
784}
785
786TH2F* Calib::get2DHistogramFromPath( TDirectoryFile* rodDir, const std::string & moduleName, const std::string & histName, int charge) {
787 std::string suffix;
788 if (charge >= 0)
789 suffix = std::format ("/C{}", charge);
790 std::string fullHistoPath = moduleName + "/" + histName + "/A0/B0" + suffix;
791 TDirectoryFile *histDir = static_cast<TDirectoryFile *>(rodDir->GetDirectory(fullHistoPath.c_str()));
792
793 if(!histDir){
794 printf("Error - Directory \"%s\" not found. Exiting..\n",fullHistoPath.c_str());
795 return nullptr;
796 }
797 TH2F *pTH2 = static_cast<TH2F*>((static_cast<TKey*>(histDir->GetListOfKeys()->First()))->ReadObj());
798 pTH2->SetDirectory(0);
799
800 return pTH2;
801}
802
803bool Calib::moduleInPart(const std::string & modName) {
804 if (modName == "DSP_ERRORS") {
805 return false;
806 }
807 return modName.starts_with(m_MODprefixes[m_whichPart]);
808}
809
810
811std::vector<float> Calib::getParams(const TF1 *f, unsigned int params){
812 std::vector<float> v;
813 for(unsigned int i = 0; i<params; i++){
814 v.push_back(f->GetParameter(i));
815 }
816
817 return v;
818}
819
820std::vector<float> Calib::getParams_quality(const TF1 *f){
821 std::vector<float> v;
822
823 v.push_back(f->GetChisquare());
824 v.push_back(f->GetNDF());
825
826 return v;
827}
828
829void Calib::graphTitles(const std::unique_ptr<TGraphErrors> &graph, const std::string &name, const std::string &Yname){
830 graph->SetTitle((std::string(name)+";Charge;"+std::string(Yname)).c_str());
831 graph->SetMarkerStyle(20);
832}
833
834
double charge(const T &p)
Definition AtlasPID.h:997
Header file for AthHistogramAlgorithm.
bool m_runOneMOD
Definition Calib.h:62
static constexpr int m_qthresh
Definition Calib.h:99
bool totFitting(const pix::PixelMapping &pm, const std::string &inTimFile, std::map< unsigned int, std::vector< std::unique_ptr< CalibFrontEndInfo > > > &map_info)
Definition Calib.cxx:16
bool m_savefile
Definition Calib.h:59
static constexpr std::array< float, m_ncharge > m_chargeArr
Definition Calib.h:95
static constexpr int m_totnbins
Definition Calib.h:85
std::vector< float > getParams(const TF1 *f, unsigned int params)
Definition Calib.cxx:811
bool fillTiming(const pix::PixelMapping &pm, const std::string &inTimFile, std::map< unsigned int, std::vector< std::unique_ptr< CalibFrontEndInfo > > > &map_info)
Definition Calib.cxx:371
TH2F * get2DHistogramFromPath(TDirectoryFile *rodDir, const std::string &moduleName, const std::string &histName, int charge=-1)
Definition Calib.cxx:786
std::vector< float > getParams_quality(const TF1 *f)
Definition Calib.cxx:820
static constexpr float m_chi_error
Definition Calib.h:65
bool reFit_normalPix(std::vector< float > &params, std::vector< float > &q, std::vector< float > &qerr, std::vector< float > &tot, std::vector< float > &toterr, std::vector< float > &sig, std::vector< float > &sigerr, const unsigned int fe)
Definition Calib.cxx:317
static constexpr float m_timHi
Definition Calib.h:83
static constexpr float m_totLo
Definition Calib.h:86
static constexpr float m_totsigHi
Definition Calib.h:90
TIter getRodIterator(const TFile &inputFile)
Definition Calib.cxx:775
static constexpr int m_phiBins
Definition Calib.h:73
static constexpr float m_sigLo
Definition Calib.h:78
static constexpr int m_ncharge
Definition Calib.h:94
static constexpr float m_totsigLo
Definition Calib.h:89
int chipId(int iphi, int ieta)
Definition Calib.cxx:738
static constexpr int m_timnbins
Definition Calib.h:81
const std::array< std::string, 4 > m_MODprefixes
Definition Calib.h:69
static constexpr int m_etaBins
Definition Calib.h:72
int pixelType(int iphi, int ieta, bool isForTOT=false)
Definition Calib.cxx:754
static constexpr float m_thrLo
Definition Calib.h:76
static constexpr int m_totsigNBins
Definition Calib.h:88
bool fillThresholds(const pix::PixelMapping &pm, const std::string &inThrFile, std::map< unsigned int, std::vector< std::unique_ptr< CalibFrontEndInfo > > > &map_info)
Definition Calib.cxx:552
int m_whichPart
Definition Calib.h:68
TIter getModuleIterator(TDirectoryFile *rodDir)
Definition Calib.cxx:781
static constexpr float m_sigHi
Definition Calib.h:79
static constexpr int m_nFE
Definition Calib.h:93
std::unique_ptr< TFile > m_wFile
Definition Calib.h:60
static constexpr std::array< float, m_ncharge > m_chargeErrArr
Definition Calib.h:96
void graphTitles(const std::unique_ptr< TGraphErrors > &graph, const std::string &name, const std::string &Yname)
Definition Calib.cxx:829
bool moduleInPart(const std::string &modName)
Definition Calib.cxx:803
static constexpr float m_totHi
Definition Calib.h:87
std::string m_testMOD
Definition Calib.h:63
static constexpr float m_timLo
Definition Calib.h:82
static constexpr int m_thrnbins
Definition Calib.h:75
static constexpr float m_thrHi
Definition Calib.h:77
bool contains(const std::string &geographicalID) const
int getID(const std::string &geographicalID) const