ATLAS Offline Software
TimingClass.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
10 #include "LArSamplesMon/Data.h"
11 #include "LArSamplesMon/OFC.h"
12 #include "LArCafJobs/Geometry.h"
13 
14 #include "TFile.h"
15 
16 using namespace TMath;
17 using namespace std;
18 
19 // *************************************** //
21 // *************************************** //
22 {
23 
24 }
25 
26 
27 
28 // **************************************** //
30 // **************************************** //
31 {
32 
33 }
34 
35 
36 // ***************************************************************************** //
37 void LArSamples::TimingClass::timePerFebAllFebs(const std::string& nrun, const std::string& name)
38 // ***************************************************************************** //
39 {
40  ofstream file;
41  string fname = "TimingFile" + nrun + "_" + name + ".txt";
42 
43  if(FileEmptyCheck(fname))
44  {
45  std::cout<<" +++++ The file:" << fname <<" is empty, time-fit will not be performed +++"<<std::endl;
46  return;
47  }
48  else
49  file.open(fname.c_str(), ios::out);
50 
51  for( uint i = 0; i < m_interface->nChannels(); i++ ){
52  if( i%10000 == 0 ) cout << "Processing entry " << i << endl;
53  if( m_interface->historySize(i) == 0 ) continue;
54  const History* history = m_interface->cellHistory(i); // extract the history of the cells of each channel
55  if( !history )continue;
56 
57  int calo = Abs(history->cellInfo()->calo());
58  int caloPart = history->cellInfo()->calo();;
59 
60  int det = 9;
61  if( caloPart > 0 ) det = 1; //sideA
62  if( caloPart < 0 ) det = 0; //sideC
63 
64  int layer = history->cellInfo()->layer();
65  int ft = history->cellInfo()->feedThrough() ;
66  int slot = history->cellInfo()->slot() ;
67 
68  for( uint j = 0; j < history->nData(); j++ ){
69  //int gain = history->data(j)->gain();
70  //int NumberRun = history->data(j)->run();
71  //remove bad channels
72  if( history->data(j)->problems(true) != "None" ) continue;
73  //noise cut
74  if( (history->data(j)->energy()/history->data(j)->noise()) < 5 ) continue;
75  //Time cut
76  if( fabs( history->data(j)->ofcTime() ) > 20.00 )continue;
77  //Energy Threshold
78  double t = history->data(j)->ofcTime(); //time
79  double e = history->data(j)->energy(); //energy
80  if( !EnergyThreshold( calo, layer, history->data(j)->quality(), ft, slot, e, t ) )continue;
81  //Fill the timing distributions with energy weighting
82  //if( calo <= 3 && gain == 0 ) h[det][ft][slot]->Fill( t,e );
83  //if( ( calo == 4 && gain == 1 ) || ( calo == 5 && gain == 1 && NumberRun < 216867 ) || ( calo == 5 && gain == 0 && NumberRun >= 216867) )h[det][ft][slot]->Fill( t,e );
84  //Print to file
85  int detec = -1;
86  if( name == "EMB" ) detec = 0;
87  else if( name == "EMEC" ) detec = 1;
88  else if( name == "HEC" ) detec = 2;
89  else if( name == "FCAL" ) detec = 3;
90  file << detec << " " << det << " " << ft << " " << slot << " " << t << " "<< e << endl;
91  }
92  }
93 
94 
95  return;
96 
97 }
98 
99 // ********************************************************************************** //
101 // ********************************************************************************** //
102 {
103  std::ifstream file(fname);
104  return file.peek() == std::ifstream::traits_type::eof();
105 }
106 
107 // ********************************************************************************** //
108 void LArSamples::TimingClass::fitTimePerFebAllFebs(const std::string& nrun, const std::string& name)
109 // ********************************************************************************** //
110 {
111 
112  for( int d = 0; d < 2; d++ ){ //side
113  for( int ft = 0; ft < 32; ft++ ){ //feedthrough
114  for( int sl = 0; sl < 15; sl++ ) //slot
115  Median[d][ft][sl] = -99.;
116  }
117  }
118 
119  TH1F *timefeb[2][32][16];
120  TF1 *fit1 = new TF1( "fit1" , "gaus" , 3 );
121 
122  double Hist_entries, Fit_mean, Fit_sigma, Fit_error;
123 
124  std::vector< std::vector<double> > myvec;
125  string tfilename;
126  tfilename = "TimingFile" + nrun + "_" + name + ".txt";
127 
128  if(FileEmptyCheck(tfilename))
129  {
130  std::cout<<" +++++ The file:" << tfilename <<" is empty, time-fit will not be performed +++"<<std::endl;
131  return ;
132  }
133 
134  myvec = readTimingFiles(tfilename);
135 
136  string Filename = "OFCTime_PerFEB_" + name + ".root";
137  TFile *FebTime = new TFile( Filename.c_str(), "READ");
138 
139  double median2, TotalE = 0.0;
140  int dete = -1, side, feedT, slot;
141 
142  std::vector<double> myvector;
143  std::vector<double> myvector2;
144  std::vector<double> energy;
145 
146  ofstream ff;
147  string ffname = "mediantest_" + nrun + "_" + name + ".txt";
148  ff.open(ffname.c_str(), ios::out);
149 
150  for( uint l = 0; l < (myvec.size() -1); l++ ){
151  float mean = myvec[l][4];
152  double E = myvec[l][5];
153  myvector.push_back(mean);
154  myvector2.push_back(mean);
155  energy.push_back(E);
156  TotalE = TotalE + E;
157  //Fill the vector when the number of slot changes or when it reaches the end of the file
158  if( ( myvec[l][3] != myvec[l+1][3] ) || ( l == myvec.size()-2 ) ){
159  //Get the median of the time for each feb (defined by its slot and ft)
160  median2 = getTimeWeightedMedian(myvector, myvector2, energy, TotalE);
161  myvector.clear();
162  myvector2.clear();
163  energy.clear();
164  TotalE = 0;
165  side = int(myvec[l][1]);
166  feedT = int(myvec[l][2]);
167  slot = int(myvec[l][3]);
168  if( name == "EMB" ) dete = 0;
169  else if( name == "EMEC" ) dete = 1;
170  else if( name == "HEC" ) dete = 2;
171  else if( name == "FCAL" ) dete = 3;
172  ff << dete << " " << side << " " << feedT << " " << slot << " " << median2 << endl;
173  Median[side][feedT][slot] = median2;
174  }
175  }
176 
177  string savename = "FEB_time_fitMean_" + nrun + "_" + name + ".txt";
178  ofstream file;
179  file.open(savename.c_str(), ios::out);
180 
181  double median_error;
182  //double mean_error;
183  std::vector<double> timevector;
184 
185  int det = -1;
186  if( name == "EMB" ) det = 0;
187  else if( name == "EMEC" ) det = 1;
188  else if( name == "HEC" ) det = 2;
189  else if( name == "FCAL" ) det = 3;
190 
191  for( int d = 0; d < 2; d++ ){ //sides
192  for( int ft = 0; ft < 32; ft++ ){ //feedthrough
193  for( int sl = 0; sl < 15; sl++ ){ //slot
194 
195  ostringstream histname;
196  histname << "h_" << name << "_" << d << "_" << ft << "_" << sl+1;
197  timefeb[d][ft][sl+1] = (TH1F*)FebTime->Get(histname.str().c_str());
198  if( !timefeb[d][ft][sl+1] )
199  break;
200 
201  Hist_entries = timefeb[d][ft][sl+1]->GetEntries();
202  //Fit only is histograms are full
203  if( Hist_entries == 0 ) continue;
204 
205  timefeb[d][ft][sl+1]->Fit( "fit1","mE","N",-10, 10 );
206 
207  //double M = fit1->GetParameter(1);
208  double rms, rmsDist;
209  rms = fit1->GetParameter(2);
210  double n1 = -99., n2 = -99.;
211  rmsDist=timefeb[d][ft][sl+1]->GetRMS();
212  double max;
213  max = timefeb[d][ft][sl+1]->GetXaxis()->GetBinCenter(timefeb[d][ft][sl+1]->GetMaximumBin());
214 
215  if( name == "FCAL" ) n1 = n2 = 0.8;
216  else if( name == "EMEC" ){
217  if( sl+1 == 13 && rmsDist > 3.0 ){
218  n1 = 1.5; n2 = 1.;
219  }
220  else if( rmsDist > 3. && sl+1 > 13 ){
221  n1 = 1.5; n2 = 0.3;
222  }
223  else if( sl+1 == 9 && ( ft == 2 || ft == 9 || ft == 15 || ft == 21 ) )
224  n1 = n2 = 3;
225  else if( ( sl+1 == 8 || sl+1 == 9 ) && rmsDist > 3.0 )
226  n1 = n2 = 2;
227  else if( sl+1 == 8 && ft == 13 )
228  n1 = n2 = 3;
229  else
230  n1 = n2 = 1.8;
231  }
232  else if( name == "HEC" ) n1 = n2 = 1.5;
233  if( name == "EMB"){
234  if( sl+1 == 9 || sl+1 == 10 ){
235  n1 = 3; n2 = 2;
236  }
237  else if( sl+1 >= 2 && sl+1 < 6 ){
238  n1 = 1.5; n2 = 2.5;
239  }
240  else if( sl+1 == 14 ){
241  n1 = 2; n2 = 1;
242  }
243  else
244  n1 = n2 = 1.8;
245  }
246 
247  timefeb[d][ft][sl+1]->Fit("fit1","mE","N", max-n1*rms, max+n2*rms);
248  Fit_mean = fit1->GetParameter(1);
249  Fit_sigma = fit1->GetParameter(2);
250  Fit_error = fit1->GetParError(1);
251 
252  median_error = 1.253*timefeb[d][ft][sl+1]->GetRMS()/sqrt(Hist_entries);
253  //mean_error = timefeb[d][ft][sl+1]->GetRMS()/sqrt(Hist_entries);
254 
255  if( Hist_entries <= 50 ){
256  file << det << " " << d << " " << ft << " " << sl+1 << " " << Median[d][ft][sl+1] << " " << median_error << endl;
257  }
258  else if( gMinuit->fCstatu == "SUCCESSFUL" && Fit_error >= Fit_mean && Fit_sigma > 1.5*rmsDist ){
259  file << det << " " << d << " " << ft << " " << sl+1 << " " << Median[d][ft][sl+1] << " " << median_error << endl;
260  }
261  else if( gMinuit->fCstatu != "SUCCESSFUL" ){
262  file << det << " " << d << " " << ft << " " << sl+1 << " " << Median[d][ft][sl+1] << " " << median_error << endl;
263  }
264  else{
265  file << det << " " << d << " " << ft << " " << sl+1 << " " << Fit_mean << " " << Fit_error << endl;
266  }
267 
268  }
269  }
270  }
271 
272 
273  ff.close();
274  file.close();
275 
276  return;
277 
278 }
279 
280 
281 /************************************************************/
282 void LArSamples::TimingClass::Time(int dete, const std::string& nrun)
283 /************************************************************/
284 {
285 
286  TH1F *h = new TH1F( Form("h_%d", dete) , Form("h_%d", dete) , 160, -20, 20 );
287  h->Sumw2();
288 
289  TH1F *h1 = new TH1F( Form("h1_%d", dete) , Form("h1_%d", dete) , 100, 0, 15000 );
290  TH2F *h2 = new TH2F( Form("h2_%d", dete) , Form("h2_%d", dete) , 100, 0, 15000, 160, -20, 20 );
291 
292 
293  for( uint i = 0; i < m_interface->nChannels(); i++ ){
294  if( m_interface->historySize(i) == 0 )continue;
295  //Extract the history of the cells of each channel
296  const History* history = m_interface->cellHistory(i);
297  if( !history )continue;
298 
299  int calo = Abs( history->cellInfo()->calo() );
300  //int caloPart = history->cellInfo()->calo();
301  int layer = history->cellInfo()->layer();
302  int slot = history->cellInfo()->slot();
303  int ft = history->cellInfo()->feedThrough();
304  //int detIndex = 1;
305  //if( calo == 1 )detIndex = 0;
306  //int side = 1;
307  //if( caloPart < 0 )side = 0;
308  //int channel = history->cellInfo()->channel();
309 
310 
311  for( uint j = 0; j < history->nData(); j++ ){
312  int Gain = history->data(j)->gain();
313  int NumberRun = history->data(j)->run();
314  //Bad Channels
315  if( history->data(j)->problems(true) != "None" ) continue;
316  //Noise cut
317  if( ( history->data(j)->energy() / history->data(j)->noise() ) < 5 )continue;
318  //Time cut
319  if( fabs( history->data(j)->ofcTime() ) > 20 )continue;
320  //Energy Threshold
321  double t = history->data(j)->ofcTime(); //time
322  double e = history->data(j)->energy(); //energy
323  if( !EnergyThreshold( calo, layer, history->data(j)->quality(), ft, slot, e, t ) )continue;
324  //Fill distributions
325  if( ( calo <= 3 && Gain == 0 ) || ( calo == 4 && Gain == 1 ) || ( calo == 5 && Gain == 1 && NumberRun < 216867 ) || ( calo == 5 && Gain == 0 && NumberRun >= 216867 ) ) h->Fill( t, e );
326  if( ( calo <= 3 && Gain == 0 ) || ( calo == 4 && Gain == 1 ) || ( calo == 5 && Gain == 1 && NumberRun < 216867 ) || ( calo == 5 && Gain == 0 && NumberRun >=216867 ) ) h1->Fill( history->data(j)->quality() );
327  if( ( calo <= 3 && Gain == 0 ) || ( calo == 4 && Gain == 1 ) || ( calo == 5 && Gain == 1 && NumberRun < 216867 ) || ( calo == 5 && Gain == 0 && NumberRun >= 216867) ) h2->Fill( history->data(j)->quality(), t );
328  }
329  }
330 
331  TCanvas *c = new TCanvas( Form("c_%d_%s", dete, nrun.c_str()), Form("c_%d_%s", dete, nrun.c_str()), 129,165,700,600 );
332  c->cd();
333  h->SetTitle(Form("hist_%d", dete) );
334  h->SetLineColor(1);
335  h->SetLineWidth(2);
336  h->GetXaxis()->SetTitle("ofcTime [ns]");
337  h->SetMarkerSize(0);
338 
339  TF1 *fit = new TF1("fit", "gaus", 3);
340  fit->SetParameter(1, h->GetMean());
341  fit->SetParameter(2, h->GetRMS());
342  fit->SetParLimits(1, -10, 10);
343  fit->SetParLimits(2, 0, 10);
344 
345  double mean, rms;
346  mean = h->GetMean();
347  rms = h->GetRMS();
348  h->Fit("fit", "", "", mean-rms, mean+0.5*rms);
349  h->Draw("histsame""E2");
350 
351  TFile *fi = new TFile( Form("Time_%d_%s.root", dete, nrun.c_str()), "RECREATE" );
352  h->Write();
353  fit->Write();
354  h1->Write();
355  h2->Write();
356  fi->Close();
357 
358  TString path = "Plots/";
359  c->Print( path+Form("Time_%d_%s.eps", dete, nrun.c_str()) );
360  c->Print( path+Form("Time_%d_%s.png", dete, nrun.c_str()) );
361 
362  return;
363 
364 }
365 
366 // ******************************************************************************* //
367 void LArSamples::TimingClass::PlotFebAverageTime24(const std::string& nrun, const std::string& name)
368 // ******************************************************************************* //
369 {
370  std::string Filename = "FEB_time_fitMean_" + nrun + "_" + name + ".txt";
371 
372  std::ifstream f(Filename.c_str(), ios::in);
373  std::vector<double> mean;
374  std::vector<int> side;
375 
376  while (!f.eof())
377  {
378  int n1, n2, n3, n4;
379  double n5, n6;
380 
381  f >> n1 >> n2 >> n3 >> n4 >> n5 >> n6;
382 
383  if (f.good() && !f.bad() && !f.fail())
384  {
385  mean.push_back(n5);
386  side.push_back(n2);
387  }
388  }
389 
390  f.close();
391 
392  TH1F *h = new TH1F(Form("FEB_Av_%s", name.c_str()), Form("FEB_Av_%s", name.c_str()), 240, -30., 30.);
393  TH1F *sAh = new TH1F(Form("FEB_Av_%s_sideA", name.c_str()), Form("FEB_Av_%s_sideA", name.c_str()), 240, -30., 30.);
394  TH1F *sCh = new TH1F(Form("FEB_Av_%s_sideC", name.c_str()), Form("FEB_Av_%s_sideC", name.c_str()), 240, -30., 30.);
395 
396  for (unsigned int i = 0; i < mean.size(); i++)
397  {
398  h->Fill(mean[i]);
399  // Fill sides
400  if (side.at(i) == 1) // A
401  sAh->Fill(mean[i]);
402  else // C
403  sCh->Fill(mean[i]);
404  }
405 
406  h->Sumw2();
407  sAh->Sumw2();
408  sCh->Sumw2();
409 
410  TCanvas *c = new TCanvas(Form("c_FEB_Av_%s", name.c_str()), Form("c_FEB_Av_%s", name.c_str()), 129, 165, 700, 600);
411  c->cd();
412  c->SetLogy();
413  h->GetXaxis()->SetTitle("<t_{FEB}> [ns]");
414  h->GetYaxis()->SetTitle("Number of FEBs / 0.25 ns");
415  h->Draw("hist");
416 
417  TString path = "Plots/";
418  c->Print(path + Form("t_FEB_%s.png", name.c_str()));
419  c->Print(path + Form("t_FEB_%s.eps", name.c_str()));
420  c->Print(path + Form("t_FEB_%s.pdf", name.c_str()));
421 
422  TCanvas *sAc = new TCanvas(Form("c_FEB_Av_%s_sideA", name.c_str()), Form("c_FEB_Av_%s_sideA", name.c_str()), 129, 165, 700, 600);
423  sAc->cd();
424  sAc->SetLogy();
425  sAh->GetXaxis()->SetTitle("<t_{FEB}> [ns]");
426  sAh->GetYaxis()->SetTitle("Number of FEBs / 0.25 ns");
427  sAh->Draw("hist");
428  sAc->Print(path + Form("t_FEB_%s_sideA.png", name.c_str()));
429 
430  TCanvas *sCc = new TCanvas(Form("c_FEB_Av_%s_sideC", name.c_str()), Form("c_FEB_Av_%s_sideC", name.c_str()), 129, 165, 700, 600);
431  sCc->cd();
432  sCc->SetLogy();
433  sCh->GetXaxis()->SetTitle("<t_{FEB}> [ns]");
434  sCh->GetYaxis()->SetTitle("Number of FEBs / 0.25 ns");
435  sCh->Draw("hist");
436  sCc->Print(path + Form("t_FEB_%s_sideC.png", name.c_str()));
437 
438  std::string plotname = "FEB_average_" + nrun + "_" + name + ".root";
439  TFile *fi = new TFile(plotname.c_str(), "RECREATE");
440  h->Write();
441  fi->Close();
442 
443  std::string plotnamesA = "FEB_average_" + nrun + "_" + name + "_sideA.root";
444  TFile *fisA = new TFile(plotnamesA.c_str(), "RECREATE");
445  sAh->Write();
446  fisA->Close();
447 
448  std::string plotnamesC = "FEB_average_" + nrun + "_" + name + "_sideC.root";
449  TFile *fisC = new TFile(plotnamesC.c_str(), "RECREATE");
450  sCh->Write();
451  fisC->Close();
452 
453  return;
454 }
455 
456 // ******************************************************************************* //
457 void LArSamples::TimingClass::PlotFebAverageTime(const std::string& nrun, const std::string& name)
458 // ******************************************************************************* //
459 {
460  string Filename = "FEB_time_fitMean_" + nrun + "_" + name + ".txt";
461 
462  ifstream f( Filename.c_str(), ios::in );
463  vector<double> mean;
464  vector<int> side;
465 
466  while( !f.eof() ){
467  int n1, n2, n3, n4;
468  double n5, n6;
469 
470  f >> n1 >> n2 >> n3 >> n4 >> n5 >> n6;
471 
472  if( f.good() && !f.bad() && !f.fail() ){
473  mean.push_back(n5);
474  side.push_back( n2 );
475  }
476  }
477 
478  f.close();
479 
480  TH1F *h = new TH1F( Form("FEB_Av_%s", name.c_str()), Form("FEB_Av_%s", name.c_str()), 120, -15., 15. );
481  TH1F *sAh = new TH1F( Form("FEB_Av_%s_sideA", name.c_str()), Form("FEB_Av_%s_sideA", name.c_str()), 120, -15., 15. );
482  TH1F *sCh = new TH1F( Form("FEB_Av_%s_sideC", name.c_str()), Form("FEB_Av_%s_sideC", name.c_str()), 120, -15., 15. );
483 
484  for( unsigned int i=0; i<mean.size(); i++ ){
485  h->Fill( mean[i] );
486  //Fill sides
487  if( side.at(i) == 1 ) //A
488  sAh->Fill( mean[i] );
489  else //C
490  sCh->Fill( mean[i] );
491  }
492 
493  h->Sumw2();
494  sAh->Sumw2();
495  sCh->Sumw2();
496 
497  TCanvas *c = new TCanvas( Form("c_FEB_Av_%s", name.c_str()), Form("c_FEB_Av_%s", name.c_str()), 129,165,700,600 );
498  c->cd();
499  c->SetLogy();
500  h->GetXaxis()->SetTitle("<t_{FEB}> [ns]");
501  h->GetYaxis()->SetTitle("Number of FEBs / 0.25 ns");
502  h->Draw("hist");
503 
504  TString path = "Plots/";
505  c->Print( path+ Form("t_FEB_%s.png", name.c_str()) );
506  c->Print( path+ Form("t_FEB_%s.eps", name.c_str()) );
507  c->Print( path+ Form("t_FEB_%s.pdf", name.c_str()) );
508 
509  TCanvas *sAc = new TCanvas( Form("c_FEB_Av_%s_sideA", name.c_str()), Form("c_FEB_Av_%s_sideA", name.c_str()), 129,165,700,600 );
510  sAc->cd();
511  sAc->SetLogy();
512  sAh->GetXaxis()->SetTitle("<t_{FEB}> [ns]");
513  sAh->GetYaxis()->SetTitle("Number of FEBs / 0.25 ns");
514  sAh->Draw("hist");
515  sAc->Print( path+ Form("t_FEB_%s_sideA.png", name.c_str()) );
516 
517  TCanvas *sCc = new TCanvas( Form("c_FEB_Av_%s_sideC", name.c_str()), Form("c_FEB_Av_%s_sideC", name.c_str()), 129,165,700,600 );
518  sCc->cd();
519  sCc->SetLogy();
520  sCh->GetXaxis()->SetTitle("<t_{FEB}> [ns]");
521  sCh->GetYaxis()->SetTitle("Number of FEBs / 0.25 ns");
522  sCh->Draw("hist");
523  sCc->Print( path+ Form("t_FEB_%s_sideC.png", name.c_str()) );
524 
525 
526  std::string plotname = "FEB_average_" + nrun + "_" + name + ".root";
527  TFile *fi = new TFile( plotname.c_str(), "RECREATE" );
528  h->Write();
529  fi->Close();
530 
531  return;
532 
533 }
534 
535 
536 // ******************************************************* //
537 void LArSamples::TimingClass::MergeFebTime( const std::string& nrun )
538 // ******************************************************* //
539 {
540  ofstream mergedfile;
541  string name = "FEB_time_fitMean_" + nrun + ".txt";
542  mergedfile.open( name.c_str(), ios::out );
543 
544  string detparts[4] = {"EMB", "EMEC", "HEC", "FCAL"};
545 
546  for( int i = 0; i < 4; i++ ){
547  string tmpname = detparts[i];
548  string file = "FEB_time_fitMean_" + nrun + "_" + tmpname + ".txt";
549 
550  if(FileEmptyCheck(file))
551  {
552  std::cout<<" +++++ no information for " << tmpname <<", therefore it's timing information is not merged +++"<<std::endl;
553  continue;
554  }
555 
556  ifstream f( file.c_str(), ios::in );
557  while( !f.eof() ){
558 
559  int n1, n2, n3, n4;
560  double n5, n6;
561 
562  f >> n1 >> n2 >> n3 >> n4 >> n5 >> n6;
563 
564  if( f.good() && !f.bad() && !f.fail() )
565  mergedfile << n1 << " " << n2 << " " << n3 << " " << n4 << " " << n5 << " " << n6 << endl;
566  }
567 
568  f.close();
569  }
570 
571  mergedfile.close();
572 
573  return;
574 
575 }
576 
577 
578 // ****************************************************************************** //
579 void LArSamples::TimingClass::getFebCorrection( const std::string& nrun )
580 // ****************************************************************************** //
581 {
582  // The structure of the Feb correction file is different than what we need
583  // The corrections are needed for each slot and not for each feed through
584  // This function rearranges the previously obtained feb correction file to put it in the right format
585 
586  for( int i = 0; i < 4; i++ ){
587  for( int j = 0; j < 2; j++ ){
588  for( int k = 0; k < 16; k++ ){
589  for( int m = 0; m < 32; m++ ){
590  param[i][j][m][k] = -99.;
591  error[i][j][m][k] = -99;
592  }
593  }
594  }
595  }
596 
597  ofstream febCorrfile;
598  string name = "FEB_Corr_" + nrun + ".txt";
599  febCorrfile.open( name.c_str(), ios::out );
600 
601  string name2 = "FEB_time_fitMean_"+nrun+".txt";
602  ifstream f( name2.c_str(), ios::in );
603 
604  if( f ){
605  while( !f.eof() ){
606  int n1, n2, n3, n4;
607  double n5, n6;
608 
609  f >> n1 >> n2 >> n3 >> n4 >> n5 >> n6;
610  if (n1 < 0 || n1 >= 4 ||
611  n2 < 0 || n2 >= 2 ||
612  n3 < 0 || n3 >= 32 ||
613  n4 < 0 || n4 >= 16)
614  {
615  cerr << "Bad line from timing file " << n1 << " " << n2 << " " << n3 << " " << n4 << endl;
616  continue;
617  }
618 
619  if( f.good() && !f.bad() && !f.fail() ){
620  param[n1][n2][n3][n4] = n5;
621  error[n1][n2][n3][n4] = n6;
622  }
623  }
624  }
625 
626  f.close();
627 
628  string HECFEBID[48] = {"0x3a1a0000", "0x3a520000", "0x3a820000", "0x3ab20000", "0x3a1a8000", "0x3a528000", "0x3a828000", "0x3ab28000", "0x3a1b0000", "0x3a530000", "0x3a830000", "0x3ab30000", "0x3a1b8000", "0x3a538000", "0x3a838000", "0x3ab38000", "0x3a1c0000", "0x3a540000", "0x3a840000", "0x3ab40000", "0x3a1c8000", "0x3a548000", "0x3a848000", "0x3ab48000", "0x3b1a0000", "0x3b520000", "0x3b820000", "0x3bb20000", "0x3b1a8000", "0x3b528000", "0x3b828000", "0x3bb28000", "0x3b1b0000", "0x3b530000", "0x3b830000", "0x3bb30000", "0x3b1b8000", "0x3b538000", "0x3b838000", "0x3bb38000", "0x3b1c0000", "0x3b540000", "0x3b840000", "0x3bb40000", "0x3b1c8000", "0x3b548000", "0x3b848000", "0x3bb48000" };
629  string FCALFEBID[28] = {"0x3a300000", "0x3a308000", "0x3a310000", "0x3a318000", "0x3a320000", "0x3a328000", "0x3a330000", "0x3a340000", "0x3a348000", "0x3a350000", "0x3a358000", "0x3a360000", "0x3a368000", "0x3a370000", "0x3b300000", "0x3b308000", "0x3b310000", "0x3b318000", "0x3b320000", "0x3b328000", "0x3b330000", "0x3b340000", "0x3b348000", "0x3b350000", "0x3b358000", "0x3b360000", "0x3b368000", "0x3b370000" };
630 
631  ofstream corrfile;
632  string nameid = "FEB_Corr_FEBID_" + nrun + ".txt";
633  corrfile.open( nameid.c_str(), ios::out );
634 
635  int kk = 0;
636  for( int i = 0; i < 4; i++ ){
637  kk = 0;
638  for( int j = 0; j < 2; j++ ){
639  for( int k = 0; k < 16; k++ ){
640  for( int m = 0; m < 32; m++ ){
641  if( param[i][j][m][k] != -99. && error[i][j][m][k] != -99.){
642  febCorrfile << i << " " << j << " " << k << " " << m << " " << param[i][j][m][k] << " " << error[i][j][m][k] << endl;
643  if( i < 2 )corrfile << "0x" << std::hex << (0x38000000|i<<25|j<<24|m<<19|(k-1)<<15) << " " << std::dec << param[i][j][m][k] << endl;
644  else if( i == 2 )corrfile << HECFEBID[kk] << " " << param[i][j][m][k] << endl;
645  else corrfile << FCALFEBID[kk] << " " << param[i][j][m][k] << endl;
646  kk++;
647  }
648  }
649  }
650  }
651  }
652 
653  febCorrfile.close();
654  corrfile.close();
655 
656  return;
657 
658 }
659 
660 
661 // *************************************** //
663 // *************************************** //
664 {
665  //Plot the average feb time per feed through and per slot
666  TGraphErrors **g = new TGraphErrors*[15];
667  TGraphErrors **g1 = new TGraphErrors*[15];
668  TGraphErrors **g2 = new TGraphErrors*[15];
669  TGraphErrors **g3 = new TGraphErrors*[15];
670  TGraphErrors **g4 = new TGraphErrors*[15];
671  TGraphErrors **g5 = new TGraphErrors*[15];
672  TGraphErrors **g6 = new TGraphErrors*[15];
673 
674  TGraph** gzero = new TGraph*[15];
675  TCanvas **c = new TCanvas*[4];
676  TLatex *l = new TLatex();
677  l->SetTextSize(0.065);
678 
679  int count = 0, count1 = 0, count2 = 0, count3 = 0, count4 = 0, count5 = 0, count6 = 0;
680  TString path = "Plots/";
681 
682  for( int det = 0; det < 4; det++ ){
683  //Plot EMEC special crates
684  int special = 1;
685  if( det == 1 )special = 2;
686  for( int crates = 0; crates < special; crates++ ){
687  string cratename = "Standard";
688  if( crates == 1 ) cratename = "Special";
689  //
690  for( int side = 0; side < 2; side++ ){
691  c[side] = new TCanvas( Form("c%d%d%s", det, side, cratename.c_str() ), Form("c%d%d%s", det,side, cratename.c_str()), 1100, 600 );
692  c[side]->Divide( 5, 3 );
693 
694  for( int i = 0; i < 15; i++ ){
695  g[i] = new TGraphErrors();
696  g1[i] = new TGraphErrors();
697  g2[i] = new TGraphErrors();
698  g3[i] = new TGraphErrors();
699  g4[i] = new TGraphErrors();
700  g5[i] = new TGraphErrors();
701  g6[i] = new TGraphErrors();
702 
703  gzero[i] = new TGraph();
704 
705  for( int j = 0; j < 32; j++ ){
706  if( det == 1 && crates == 0 && ( j == 2 || j == 3 || j == 9 || j == 10 || j == 15 || j == 16 || j ==21 || j == 22 ) )continue;
707  if( det == 1 && crates == 1 && !( j == 2 || j == 3 || j == 9 || j == 10 || j == 15 || j == 16 || j ==21 || j == 22 ) )continue;
708 
709  gzero[i]->SetPoint( j, j, 0 );
710  if( param[det][side][j][i+1] == -99 )continue;
711  //
712  count++;
713  g[i]->SetPoint( count-1, j, param[det][side][j][i+1] );
714  g[i]->SetPointError( count-1, 0, error[det][side][j][i+1] );
715  //
716  if( det == 0 ){
717  g[i]->SetMarkerColor(4);
718  g[i]->GetYaxis()->SetRangeUser(-5, 5);
719  }
720  else{
721  if( j == 6 ){
722  count1++;
723  g1[i]->SetPoint(count1-1, j, param[det][side][j][i+1]);
724  g1[i]->SetPointError(count1-1, 0, error[det][side][j][i+1]);
725  g1[i]->SetMarkerColor(4);
726  }
727  else if( j == 3 || j == 10 || j == 16 || j == 22 ){
728  if( i >= 4 && i <= 9 ){
729  count2++;
730  g2[i]->SetPoint(count2-1, j, param[det][side][j][i+1]);
731  g2[i]->SetPointError(count2-1, 0, error[det][side][j][i+1]);
732  g2[i]->SetMarkerColor(3);
733  }
734  else if( i == 0 || i == 1){
735  count3++;
736  g3[i]->SetPoint(count3-1, j, param[det][side][j][i+1]);
737  g3[i]->SetPointError(count3-1, 0, error[det][side][j][i+1]);
738  g3[i]->SetMarkerColor(6);
739  }
740  }
741  else if( j == 2 || j == 9 || j == 15 || j == 21 ){
742  count4++;
743  g4[i]->SetPoint(count4-1, j, param[det][side][j][i+1]);
744  g4[i]->SetPointError(count4-1, 0, error[det][side][j][i+1]);
745  g4[i]->SetMarkerColor(2);
746  }
747  else if( i <= 12 ){
748  if( j%2 ){
749  count5++;
750  g5[i]->SetPoint(count5-1, j, param[det][side][j][i+1]);
751  g5[i]->SetPointError(count5-1, 0, error[det][side][j][i+1]);
752  g5[i]->SetMarkerColor(1+12*(j%2));
753  }
754  else{
755  count6++;
756  g6[i]->SetPoint(count6-1, j, param[det][side][j][i+1]);
757  g6[i]->SetPointError(count6-1, 0, error[det][side][j][i+1]);
758  g6[i]->SetMarkerColor(1+12*(j%2));
759  }
760  }
761  g[i]->GetYaxis()->SetRangeUser(-5, 5);
762  }
763  }
764 
765  count = 0; count1 = 0; count2 = 0; count3 = 0; count4 = 0; count5 = 0; count6 = 0;
766 
767  c[side]->cd(i+1);
768  g[i]->GetXaxis()->SetTitle("FT");
769  g[i]->GetYaxis()->SetTitle("ofcTime[ns]");
770  g[i]->GetXaxis()->SetLabelSize(0.065);
771  g[i]->GetYaxis()->SetLabelSize(0.065);
772  g[i]->SetMarkerSize(0.4);
773  g1[i]->SetMarkerSize(0.4);
774  g2[i]->SetMarkerSize(0.4);
775  g3[i]->SetMarkerSize(0.4);
776  g4[i]->SetMarkerSize(0.4);
777  g5[i]->SetMarkerSize(0.4);
778  g6[i]->SetMarkerSize(0.4);
779  g[i]->SetTitle( Form("#color[2]{#bf{Slot %d}}", i+1 ) );
780  g[i]->Draw("AP");
781  g1[i]->Draw("P");
782  g2[i]->Draw("P");
783  g3[i]->Draw("P");
784  g4[i]->Draw("P");
785  g5[i]->Draw("P");
786  g6[i]->Draw("P");
787 
788  gzero[i]->SetMarkerSize(0.2);
789  gzero[i]->Draw("P");
790  if( i == 14 ){
791  string detname = "EMB";
792  if( det == 1 )detname = "EMEC";
793  else if( det == 2 )detname = "HEC";
794  else if( det == 3 )detname = "FCAL";
795  c[side]->Print( path+ Form("t_FEB_side%d_%s_%s.png", side, detname.c_str(), cratename.c_str()));
796  }
797  }
798  }
799  }
800  }
801 
802  return;
803 
804 }
805 
806 
807 // **************************************************************************************************************************** //
808 bool LArSamples::TimingClass::EnergyThreshold( int calo, int layer, int quality, int ft, int slot, double energy, double time )
809 // **************************************************************************************************************************** //
810 {
811  bool pass = true;
812  if( quality > 4000 ) pass = false;
813 
814  //
815  if( layer == 0 )
816  if( fabs( time ) > 10.00 )pass = false;
817 
818  //
819  if( calo == 1 ){
820  if( layer == 0 && energy < 1000. ) pass = false;
821  if( layer == 1 && energy < 1000. ) pass = false;
822  if( layer == 2 && energy < 3000. ) pass = false;
823  if( layer == 3 && energy < 1500. ) pass = false;
824  }
825 
826  //
827  if( calo == 2 || calo == 3 ){
828  if( layer == 0 && energy < 1500. ) pass = false;
829  if( layer == 1 ){
830  if( calo == 3 && energy < 3000. ) pass = false;
831  if( calo == 2 && energy < 1000. ) pass = false;
832  }
833 
834  if( layer == 2 ){
835  if( calo == 3 && energy < 2000. ) pass = false;
836  if( calo == 2 && energy < 3000. ) pass = false;
837  }
838 
839  if( layer == 3 && energy < 2000. ) pass = false;
840  if( layer == 3 && quality > 10000 ) pass = false;
841  if( slot == 13 && ( ft == 2 || ft == 9 || ft == 15 || ft == 21 ) && quality > 4000 ) pass = false;
842  if( slot == 14 && ( ft == 2 || ft == 9 || ft == 15 || ft == 21 ) && quality > 4000 ) pass = false;
843  if( slot == 15 && ( ft == 2 || ft == 9 || ft == 15 || ft == 21 ) && quality > 4000 ) pass = false;
844  }
845 
846  if( calo == 4 && energy < 3500. ) pass = false;
847  if( calo == 5 && energy < 10000. ) pass = false;
848 
849  return pass;
850 
851 }
852 
853 // ******************************************************************************* //
854 vector< vector<double> > LArSamples::TimingClass::readTimingFiles(const std::string& file)
855 // ******************************************************************************* //
856 {
857  std::vector< std::vector <double> > Data;
858 
859  ifstream f( file.c_str() );
860 
861  TH1F * h[2][32][16];
862 
863  string name;
864  if( file.find( "EMB" ) < file.length() )name = "EMB";
865  else if( file.find( "EMEC" ) < file.length() )name = "EMEC";
866  else if( file.find( "HEC" ) < file.length() )name = "HEC";
867  else if( file.find( "FCAL" ) < file.length() )name = "FCAL";
868 
869  for( int d = 0; d < 2; d++ ){ //sides
870  for( int ft = 0; ft < 32; ft++ ){ //feedthrough
871  for( int sl = 0; sl < 15; sl++ ){ //slots
872  h[d][ft][sl+1] = new TH1F( Form("h_%s_%d_%d_%d", name.c_str(), d, ft, sl+1) , Form("h_%s_%d_%d_%d", name.c_str(),d,ft, sl+1), 200, -20, 20 );
873  h[d][ft][sl+1]->Sumw2();
874  }
875  }
876  }
877 
878 
879 
880  if( !f )cerr << "Timing File " << file << "cannot be read" << endl;
881  else {
882  while( !f.eof() ){
883  int n1, n2, n3, n4;
884  double n5, n6;
885 
886  if( f.good() && !f.bad() && !f.fail() ){
887  f >> n1 >> n2 >> n3 >> n4 >> n5 >> n6;
888  if (n2 < 0 || n2 >= 2 ||
889  n3 < 0 || n3 >= 32 ||
890  n4 < 0 || n4 >= 16)
891  {
892  cerr << "Bad line from timing file " << n2 << " " << n3 << " " << n4 << endl;
893  continue;
894  }
895  //
896  Data.push_back(std::vector<double>{static_cast<double>(n1),
897  static_cast<double>(n2),
898  static_cast<double>(n3),
899  static_cast<double>(n4),
900  static_cast<double>(n5),
901  static_cast<double>(n6)});
902  //
903  h[n2][n3][n4]->Fill( n5, n6 );
904  }
905  }
906  }
907 
908  f.close();
909 
910 
911  //Save the Histograms
912  string filename = "OFCTime_PerFEB_" + name + ".root";
913  TFile *ff = new TFile( filename.c_str(), "RECREATE" );
914 
915  for( int d = 0; d < 2; d++ ){ //side
916  for( int ft = 0; ft < 32; ft++ ){ //feedthrough
917  for( int sl = 0; sl < 15; sl++ ) //slot
918  h[d][ft][sl+1]->Write();
919  }
920  }
921 
922  ff->Close();
923 
924 
925 
926  return Data;
927 
928 }
929 
930 
931 // ********************************************************************************************************************************************** //
932 double LArSamples::TimingClass::getTimeWeightedMedian(std::vector<double> time, const std::vector<double>& time2, const std::vector<double>& weight, double totalW)
933 // ********************************************************************************************************************************************* //
934 {
935  TGraph *g = new TGraph();
936  double wmedian = -99., weights, w0 = -1., cumulWeights = 0.0;
937  int Size = time.size();
938  bool wcondition = false;
939 
940  sort( time.begin(), time.end() );
941 
942  for( int i = 0; i < Size; i++ ){
943  int pos = std::find( time2.begin(),time2.end(), time[i] ) - time2.begin();
944  weights = weight[pos] / totalW;
945  cumulWeights = cumulWeights + weights;
946  g->SetPoint( g->GetN(), cumulWeights, time[i] );
947  if( i == 0 ) w0 = weight[pos] / totalW;
948  wcondition = w0 > 0.5;
949  }
950 
951  if( wcondition ) wmedian = g->Eval(w0+(1-w0)/2);
952  else wmedian= g->Eval(0.5);
953 
954  g->Delete();
955 
956  return wmedian;
957 
958 }
959 
LArNewCalib_Delay_OFC_Cali.Gain
Gain
Definition: LArNewCalib_Delay_OFC_Cali.py:87
TreeShapeErrorGetter.h
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
ClusterMessageType::Data
@ Data
LArSamples::Data::noise
double noise() const
Definition: Data.cxx:64
PlotCalibFromCool.ft
ft
Definition: PlotCalibFromCool.py:329
LArSamples::TimingClass::PlotFebAverageTime24
void PlotFebAverageTime24(const std::string &nrun, const std::string &name)
Definition: TimingClass.cxx:367
AbsShapeErrorGetter.h
mean
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:254
LArSamples::CellInfo::slot
short slot() const
Definition: CellInfo.h:68
LArSamples::TimingClass::TimingClass
TimingClass()
Definition: TimingClass.cxx:20
LArSamples::CellInfo::calo
CaloId calo() const
Definition: CellInfo.h:50
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
LArSamples::History::cellInfo
const CellInfo * cellInfo() const
Definition: History.h:56
TRTCalib_Extractor.det
det
Definition: TRTCalib_Extractor.py:36
hist_file_dump.d
d
Definition: hist_file_dump.py:142
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
TrigConf::JetWindowSize::Size
Size
Definition: TriggerThresholdValue.h:17
LArSamples::History
Definition: History.h:35
LArSamples::CellInfo::feedThrough
short feedThrough() const
Definition: CellInfo.h:65
Geometry.h
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:70
LArSamples::Data::run
int run() const
Definition: Data.cxx:27
LArSamples::TimingClass::EnergyThreshold
bool EnergyThreshold(int calo, int layer, int quality, int ft, int slot, double energy, double time)
Definition: TimingClass.cxx:808
OFC.h
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:157
LArSamples::Data::ofcTime
double ofcTime() const
Definition: Data.h:111
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
read_hist_ntuple.h1
h1
Definition: read_hist_ntuple.py:21
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
LArSamples::TimingClass::MergeFebTime
void MergeFebTime(const std::string &nrun)
Definition: TimingClass.cxx:537
std::sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition: DVL_algorithms.h:554
XMLtoHeader.count
count
Definition: XMLtoHeader.py:84
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:190
TRT::Hit::side
@ side
Definition: HitInfo.h:83
LArSamples::TimingClass::fitTimePerFebAllFebs
void fitTimePerFebAllFebs(const std::string &nrun, const std::string &name)
Definition: TimingClass.cxx:108
LArSamples::TimingClass::getFebCorrection
void getFebCorrection(const std::string &nrun)
Definition: TimingClass.cxx:579
LArSamples::History::data
const Data * data(unsigned int i) const
Definition: History.cxx:91
python.changerun.kk
list kk
Definition: changerun.py:39
uint
unsigned int uint
Definition: LArOFPhaseFill.cxx:20
PixelAthClusterMonAlgCfg.histname
histname
Definition: PixelAthClusterMonAlgCfg.py:106
ParticleGun_FastCalo_ChargeFlip_Config.energy
energy
Definition: ParticleGun_FastCalo_ChargeFlip_Config.py:78
fitman.g1
g1
Definition: fitman.py:619
lumiFormat.i
int i
Definition: lumiFormat.py:85
D3PDSizeSummary.ff
ff
Definition: D3PDSizeSummary.py:305
h
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
LArSamples::Data::quality
double quality() const
Definition: Data.h:114
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
LArSamples::TimingClass::getTimeWeightedMedian
double getTimeWeightedMedian(std::vector< double > time, const std::vector< double > &time2, const std::vector< double > &weight, double totalW)
Definition: TimingClass.cxx:932
file
TFile * file
Definition: tile_monitor.h:29
mc.g4
g4
Definition: mc.PhPy8EG_A14NNPDF23_DY_VLQ_example.py:25
hist_file_dump.f
f
Definition: hist_file_dump.py:140
LArSamples::TimingClass::Time
void Time(int dete, const std::string &nrun)
Definition: TimingClass.cxx:282
LArSamples::History::nData
unsigned int nData() const
Definition: History.h:51
LArSamples::TimingClass::~TimingClass
~TimingClass()
Definition: TimingClass.cxx:29
MuonValidation_CreateResolutionProfiles.fit
def fit(h, emin, emax)
Definition: MuonValidation_CreateResolutionProfiles.py:69
fitman.g2
g2
Definition: fitman.py:624
LArSamples::Data
Definition: Data.h:72
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
LArSamples::Data::gain
CaloGain::CaloGain gain() const
Definition: Data.h:85
weights
Definition: herwig7_interface.h:38
LArSamples::TimingClass::readTimingFiles
std::vector< std::vector< double > > readTimingFiles(const std::string &file)
Definition: TimingClass.cxx:854
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
TimingClass.h
LArSamples::TimingClass::FileEmptyCheck
bool FileEmptyCheck(const std::string &fname)
Definition: TimingClass.cxx:100
Rtt_histogram.n1
n1
Definition: Rtt_histogram.py:21
python.AthDsoLogger.fname
string fname
Definition: AthDsoLogger.py:66
LArSamples::Data::problems
TString problems(bool sayNone=false) const
Definition: Data.cxx:135
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
CaloSwCorrections.time
def time(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:242
python.TrackLeptonConfig.quality
quality
Definition: TrackLeptonConfig.py:16
LArSamples::TimingClass::PlotFebtime
void PlotFebtime()
Definition: TimingClass.cxx:662
beamspotnt.rms
rms
Definition: bin/beamspotnt.py:1265
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:23
Data.h
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
History.h
LArSamples::TimingClass::PlotFebAverageTime
void PlotFebAverageTime(const std::string &nrun, const std::string &name)
Definition: TimingClass.cxx:457
error
Definition: IImpactPoint3dEstimator.h:70
Interface.h
python.compressB64.c
def c
Definition: compressB64.py:93
LArSamples::CellInfo::layer
short layer() const
Definition: CellInfo.h:53
LArSamples::TimingClass::timePerFebAllFebs
void timePerFebAllFebs(const std::string &nrun, const std::string &name)
Definition: TimingClass.cxx:37
LArSamples::Data::energy
double energy() const
Definition: Data.h:108
fitman.k
k
Definition: fitman.py:528
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106