Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
T0CalibrationClassic.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "MdtCalibT0/T0CalibrationClassic.h"
6 
7 #include <algorithm>
8 #include <iostream>
9 #include <string>
10 #include <vector>
11 #include <format>
12 
13 #include "MdtCalibT0/T0CalibrationOutput.h"
14 #include "MuonCalibEventBase/MuonCalibSegment.h"
15 #include "MuonCalibIdentifier/MuonFixedId.h"
16 #include "TH1F.h"
17 #include "TH2F.h"
18 #include "TMinuit.h"
19 
20 namespace MuonCalib {
21 
22  // function to be fitted to time spectrum
23 
24  inline Double_t TimeSpectrum_func(Double_t *xx, Double_t *par) {
25  Double_t &x(xx[0]);
26  return par[0] + (par[1] * (1 + par[2] * std::exp(-(x - par[4]) / par[3]))) /
27  ((1 + std::exp((-x + par[4]) / par[6])) * (1 + std::exp((x - par[5]) / par[7])));
28  }
29 
30  T0CalibrationClassic::T0CalibrationClassic(const std::string &name, const T0ClassicSettings *settings) :
31  IMdtCalibration(name), m_settings(settings), m_converged(false), m_name(name), m_result(nullptr), m_delete_settings(false) {
32  if (!m_settings) {
33  std::array<double, 8> params{}; // warning: this would give a memory leak
34  params[0] = 0.; // noise level
35  params[1] = 6.5;
36  params[2] = 6.5;
37  params[3] = 155.;
38  params[4] = 280.; // t0
39  params[5] = 980.; // tmax
40  params[6] = 9.; // t0 slope
41  params[7] = 11.5; // tmax slope
42  m_settings = new T0ClassicSettings(0., 300., 100, -100., 900., 1000, 1, 1000, 1, 8, params, 10., 4);
43  m_delete_settings = true;
44  }
45  MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
46  if (log.level() <= MSG::INFO) log << MSG::INFO << "T0CalibrationClassic::T0CalibrationClassic" << m_name << " " << name << endmsg;
47  if (m_settings->printLevel() <= 2) m_settings->print();
48  m_currentItnum = 1;
49  std::string HistoFileName = "T0Classic_" + m_name + ".root";
50  if (log.level() <= MSG::INFO)
51  log << MSG::INFO << "T0CalibrationClassic::T0CalibrationClassic" << m_name << " " << name << " " << HistoFileName << endmsg;
52  m_file = std::make_unique<TFile>(HistoFileName.c_str(), "recreate");
53  m_regiondir = m_file->mkdir(m_name.c_str());
54  }
55 
56  T0CalibrationClassic::~T0CalibrationClassic() {
57  m_file->Write();
58  m_file->Close();
59  if (m_delete_settings) delete m_settings;
60  }
61 
62  bool T0CalibrationClassic::handleSegment(MuonCalibSegment &seg) {
63  for (const MuonCalibSegment::MdtHitPtr& hit : seg.mdtHOT()) {
64  float distanceToRO = hit->distanceToReadout();
65 
66  bool ROside = distanceToRO < 130000.; // this means that there is no selection along the tube
67  if (ROside) {
68  const MuonFixedId &id = hit->identify();
69  const MdtIdHelper& idHelper{m_result->idHelperSvc()->mdtIdHelper()};
70  // get the T0 originally subtracted for this hit
71  int nML = id.mdtMultilayer();
72  int nL = id.mdtTubeLayer();
73  int nT = id.mdtTube();
74  const Identifier tubeId = idHelper.channelID(id.stationNameString(), id.eta(), id.phi(),
75  nML, nL, nT);
76  const MdtTubeFitContainer::SingleTubeCalib *stc = m_result->getCalib(tubeId);
77  if (!stc) {
78  MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
79  log << MSG::WARNING << "no Single Tube Calib info found for ML=" << nML << " L=" << nL << " T=" << nT << endmsg;
80  log << MSG::WARNING << "container size " << m_result->size() << endmsg;
81  log << MSG::WARNING << "container nML " << m_result->numMultilayers() << endmsg;
82  log << MSG::WARNING << "container nL " << m_result->numLayers() << endmsg;
83  log << MSG::WARNING << "container nT " << m_result->numTubes() << endmsg;
84  }
85 
86  // get histos
87  T0ClassicHistos *histos = getHistos(id.getIdInt());
88 
89  // fill histos
90 
91  float ttime = hit->driftTime();
92 
93  histos->time->Fill(ttime);
94  histos->adc->Fill(hit->adcCount());
95  // book an additional dummy set of histograms to fit the global T0
96  T0ClassicHistos *histosAll = getHistos(0);
97  histosAll->time->Fill(ttime);
98  histosAll->adc->Fill(hit->adcCount());
99 
100  // M.I. Jun20-07 ---- Adding MultiLayer Histos :
101  T0ClassicHistos *histosML = getHistos(nML);
102  histosML->time->Fill(ttime);
103  histosML->adc->Fill(hit->adcCount());
104  // M.I. Jan1507 ---- Adding Mezzanine Histos :
105  T0ClassicHistos *histosMezz = getHistos(id.mdtMezzanine());
106  histosMezz->time->Fill(ttime);
107  histosMezz->adc->Fill(hit->adcCount());
108 
109  // M.I. 16Oct07 ---- Adding "SerialGas" Histos :
110  int serialGas = nML * 10 + (nT - 1) % 3 + 1;
111  T0ClassicHistos *histosSerialGas = getHistos(serialGas);
112  histosSerialGas->time->Fill(ttime);
113  histosSerialGas->adc->Fill(hit->adcCount());
114  } // end if (ROside)
115  } // end loop over seg.mdtHOT
116  return true;
117  }
118 
119  bool T0CalibrationClassic::analyse() {
120  MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
121  if (log.level() <= MSG::INFO) log << MSG::INFO << "T0CalibrationClassic::analyse iteration " << m_currentItnum << endmsg;
122 
123  const MdtIdHelper& idHelper{m_result->idHelperSvc()->mdtIdHelper()};
124  // loop over m_histos histograms
125  for (std::unique_ptr<T0ClassicHistos> &hist : m_histos) {
126  if (m_settings->fitTime()) {
127  MdtTubeFitContainer::SingleTubeFit full;
128  MdtTubeFitContainer::SingleTubeCalib st;
129  int idtube = hist->id;
130  if (idtube != 0 && (int)(idtube / 100000000) != 9) {
131  doTimeFit(*hist, full, st);
132  doAdcFit(*hist, full, st);
133  MuonFixedId fId(idtube);
134  int nML = fId.mdtMultilayer();
135  int nL = fId.mdtTubeLayer();
136  int nT = fId.mdtTube();
137  const Identifier tubeId = idHelper.channelID(fId.stationNameString(), fId.eta(), fId.phi(),
138  nML, nL, nT);
139 
140  bool setInfo = m_result->setCalib(std::make_unique<MdtTubeFitContainer::SingleTubeCalib>(st), tubeId, log);
141  if (!setInfo) log << MSG::WARNING << "T0CalibrationClassic::PROBLEM! could not set SingleTubeCalib info " << endmsg;
142  setInfo = m_result->setFit(std::move(full), tubeId, log);
143  if (!setInfo) log << MSG::WARNING << "T0CalibrationClassic::PROBLEM! could not set SingleTubeFullInfo info " << endmsg;
144  }
145  }
146  }
147 
148  m_currentItnum++;
149  m_file->Write();
150  return true;
151  }
152 
153  T0CalibrationClassic::MdtCalibOutputPtr T0CalibrationClassic::analyseSegments(const MuonSegVec &segs) {
154  for (const auto & seg : segs) handleSegment(*seg);
155  analyse();
156  return getResults();
157  }
158 
159  void T0CalibrationClassic::doTimeFit(T0ClassicHistos &T0h, MdtTubeFitContainer::SingleTubeFit &fi,
160  MdtTubeFitContainer::SingleTubeCalib &stc) {
161  std::unique_ptr<TMinuit> gMinuit = std::make_unique<TMinuit>();
162  int fitMezz(1);
163 
164  const int np = m_settings->numParams();
165  std::vector<double> pfit(np, 0), errfit(np, 0);
166  Double_t **matrix = new Double_t *[np];
167  // initialize everything to make Coverity happy
168  for (int i = 0; i < np; i++) {
169  matrix[i] = new Double_t[np];
170  for (int j = 0; j < np; j++) matrix[i][j] = 0.;
171  }
172  const std::array<double, 8> &pdefault = m_settings->params();
173  double chi2{0};
174  int ndof{0};
175 
176  MuonFixedId fId(T0h.id);
177  int isMultilayer(0);
178 
179  if (T0h.id == 1 || T0h.id == 2 || (T0h.id > 10 && T0h.id <= 23)) isMultilayer = 1;
180 
181  MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
182  if ((fId.isValid() && fId.is_mdt()) || isMultilayer) {
183  if (log.level() <= MSG::INFO) log << MSG::INFO << " STARTING doTimeFit " << endmsg;
184  TH1 *h = nullptr;
185 
186  if (isMultilayer) {
187  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << T0h.id << endmsg;
188  h = T0h.time.get();
189  }
190 
191  if (fitMezz == 1 && !isMultilayer) { // FIT MEZZANINE
192  int hIdMezz = fId.mdtMezzanine();
193 
194  std::string HistoId = std::format("time_mezz_{:}", (hIdMezz) % (900000000));
195  if (log.level() <= MSG::DEBUG) {
196  log << MSG::DEBUG << " doTimeFit HistogramId : " << T0h.id << endmsg;
197  log << MSG::DEBUG << " doTimeFit Histogram : " << HistoId << endmsg;
198  }
199  TH1F *timeHis = (TH1F *)m_regiondir->Get(HistoId.c_str());
200  if (!timeHis) {
201  for (int i = 0; i < np; i++) delete[] matrix[i];
202  delete[] matrix;
203  return;
204  }
205 
206  T0ClassicHistos *histosMezz = getHistos(fId.mdtMezzanine());
207  h = histosMezz->time.get();
208  }
209 
210  Stat_t entries = h->GetEntries();
211  fi.statistics = (int)entries;
212 
213  if (log.level() <= MSG::VERBOSE)
214  log << MSG::VERBOSE << " histogram " << h->GetName() << " " << h->GetTitle() << " entries=" << h->GetEntries()
215  << " min entries=" << m_settings->entries() << std::endl;
216 
217  // CHECK whether the Selected Histogram has enough entries
218  if ((int)entries > m_settings->entries()) {
219  // CHECK whether the histogram has been already fitted
220  std::unique_ptr<TF1> FitFunction{h->GetFunction("TimeSpectrum")};
221  if (FitFunction) {
222  for (int i = 0; i < np; i++) {
223  pfit[i] = FitFunction->GetParameter(i);
224  errfit[i] = FitFunction->GetParError(i);
225  }
226  chi2 = FitFunction->GetChisquare(); // total chi2
227  ndof = FitFunction->GetNDF(); // number of degrees of freedom
228  } else { // The Selected Histo has NOT been fitted. Fit Now
229  std::unique_ptr<TF1> TimeSpectrum = std::make_unique<TF1>(
230  "TimeSpectrum", "[0]+([1]*(1+[2]*exp(-(x-[4])/[3])))/((1+exp((-x+[4])/[6]))*(1+exp((x-[5])/[7]))) ",
231  m_settings->minTime(), m_settings->maxTime());
232  for (int i = 0; i < np; i++) {
233  pfit[i] = pdefault[i];
234  if (log.level() <= MSG::DEBUG)
235  log << MSG::DEBUG << "T0CalibrationClassic::doTimeFit initial parameters" << i << "=" << pfit[i] << endmsg;
236  }
237  searchParams(h, &pfit[0], np);
238  if (log.level() <= MSG::DEBUG) {
239  log << MSG::DEBUG << "T0CalibrationClassic::doTimeFit parameters after searchParams " << endmsg;
240  for (int i = 0; i < np; ++i) { log << MSG::DEBUG << "i,pfit(i) " << i << " " << pfit[i] << endmsg; }
241  }
242  TimeSpectrum->SetParameters(pfit.data());
243  TimeSpectrum->SetParLimits(0, 0., 5.);
244  TimeSpectrum->SetParLimits(1, 0., 1000.);
245  TimeSpectrum->SetParLimits(2, 0., 40.);
246  TimeSpectrum->SetParLimits(3, 50., 400.);
247  TimeSpectrum->SetParLimits(4, 0., 1000.);
248  // 5 parameters fit
249  TimeSpectrum->SetParLimits(5, pfit[5], pfit[5]);
250  TimeSpectrum->SetParLimits(6, pfit[6], pfit[6]);
251  TimeSpectrum->SetParLimits(7, pfit[7], pfit[7]);
252  h->Fit("TimeSpectrum", "QLB");
253  // 6 parameters fit
254  TimeSpectrum->SetParLimits(5, 500., 2000.);
255  TimeSpectrum->SetParLimits(6, pfit[6], pfit[6]);
256  TimeSpectrum->SetParLimits(7, pfit[7], pfit[7]);
257  h->Fit("TimeSpectrum", "QLB");
258  // 7 parameters fit
259  TimeSpectrum->SetParLimits(6, 4., 30.);
260  h->Fit("TimeSpectrum", "QLB");
261  // final 8 parameters fit
262  TimeSpectrum->SetParLimits(6, 4., 30.);
263  TimeSpectrum->SetParLimits(7, 4., 30.);
264  double xmin = h->GetBinLowEdge(1);
265  double xmax = pfit[5] + 250.;
266  h->Fit("TimeSpectrum", "QLB", "", xmin, xmax);
267 
268  gMinuit->mnemat(&matrix[0][0], np);
269  for (int i = 0; i < np; i++) {
270  pfit[i] = TimeSpectrum->GetParameter(i);
271  errfit[i] = TimeSpectrum->GetParError(i);
272  }
273  chi2 = TimeSpectrum->GetChisquare(); // total chi2
274  ndof = TimeSpectrum->GetNDF(); // number of degrees of freedom
275  }
276  // THE NEW HISTOGRAM HAS BEEN FITTED
277  if (ndof == 0.) ndof = -1;
278 
279  if (log.level() <= MSG::VERBOSE)
280  log << MSG::VERBOSE << " fit results chi2/ndof=" << chi2 / ndof << " T0=" << pfit[4] << " err=" << errfit[4] << endmsg;
281  if (chi2 / ndof < m_settings->chi2max()) {
282  stc.statusCode = 0; // success
283  } else {
284  stc.statusCode = 3; // bad chi2
285  }
286  stc.t0 = pfit[4];
287  fi.chi2Tdc = chi2 / ndof;
288  for (int i = 0; i < np; i++) fi.par[i] = pfit[i];
289 
290  // NOW we get rid of the covariance matrix
291  /*******************************************************
292 
293  *******************************************************/
294  // NOW we set the the first 8 values of fi.cov to errfit
295  for (int i = 0; i < np; i++) fi.cov[i] = errfit[i];
296 
297  //
298  // Try 4 parameters fit now
299  // NOW COMMENTED BUT SHOULD BE TRIED !!!! (as it was in calib )
300  /*****************************************************
301  double lowt;
302  double hight;
303  double * pfd=new double[4];
304  pfd[0]=pfit[0];
305  pfd[1]=pfit[1];
306  pfd[2]=pfit[4];
307  pfd[3]=pfit[7];
308  lowt=pfit[4]-2.*pfit[6];
309  hight=pfit[4]+50.;
310  TF1 * FermiDirac = new TF1("FermiDirac",
311  "[0]+([1]/(1+exp(-(x-[3])/[2])))",
312  lowt,hight);
313  // FermiDirac->SetParameters(pfd);
314  //h->Fit("FermiDirac","LBV");
315  *****************************************************/
316 
317  } else {
318  stc.statusCode = 2; // too few entries
319  }
320  }
321  for (int i = 0; i < np; i++) delete[] matrix[i];
322  delete[] matrix;
323 
324  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " ENDING doTimeFit " << endmsg;
325  }
326 
327  void T0CalibrationClassic::doAdcFit(T0ClassicHistos &T0h, MdtTubeFitContainer::SingleTubeFit &fi,
328  MdtTubeFitContainer::SingleTubeCalib &stc) {
329  // THIS SETTINGS VARIABLE HAS TO BE IMPLEMENTED (as of March 7, 2007) :
330  // int fitMezz = m_settings->fitMezzanine();
331  // int fitMezz(123) ;
332  int fitMezz(1);
333 
334  double chi2(0);
335  int ndof(0);
336  std::array<double, 4> par{0}, errpar{0};
337 
338  MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
339  MuonFixedId fId(T0h.id);
340  if (fId.isValid() && fId.is_mdt()) {
341  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " doAdcFit : checking Single tube entries" << endmsg;
342  double adcThreshold = 50.;
343  TH1 *hcheck = T0h.adc.get();
344  int nhits = static_cast<int>(hcheck->GetEntries());
345 
346  int adcBinThreshold = (int)((adcThreshold - hcheck->GetBinLowEdge(1)) / (hcheck->GetBinWidth(1)) + 1); //
347  int nhitsAboveAdcCut = (int)hcheck->Integral(adcBinThreshold, hcheck->GetNbinsX());
348  if (log.level() <= MSG::DEBUG)
349  log << MSG::DEBUG << " doAdcFit : TotHits, nhitsAboveAdcCut " << nhits << " " << nhitsAboveAdcCut << endmsg;
350 
351  fi.cov[20] = nhits;
352  fi.cov[21] = nhitsAboveAdcCut;
353 
354  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " STARTING doAdcFit " << endmsg;
355 
356  TH1 *h = nullptr;
357 
358  if (fitMezz == 1) {
359  int hIdMezz = fId.mdtMezzanine();
360  std::string HistoId = std::format("charge_mezz_{:}", (hIdMezz) % (900000000));
361  if (log.level() <= MSG::DEBUG) {
362  log << MSG::DEBUG << " doAdcFit HistogramId : " << T0h.id << endmsg;
363  log << MSG::DEBUG << " doAdcFit Histogram : " << HistoId << endmsg;
364  }
365  TH1F *adcHis = (TH1F *)m_regiondir->Get(HistoId.c_str());
366  if (!adcHis) return;
367 
368  T0ClassicHistos *histosMezz = getHistos(fId.mdtMezzanine());
369  h = histosMezz->adc.get();
370  }
371 
372  if (log.level() <= MSG::VERBOSE)
373  log << MSG::VERBOSE << " histogram " << h->GetName() << " " << h->GetTitle() << " entries=" << h->GetEntries() << endmsg;
374  Stat_t entries = h->GetEntries();
375 
376  // CHECK whether the Selected Histogram has enough entries
377  if ((int)entries > m_settings->entries()) {
378  // CHECK whether the histogram has been already fitted
379  TF1 *FitFunction = h->GetFunction("AdcSpectrum");
380  if (FitFunction) {
381  chi2 = FitFunction->GetChisquare();
382  ndof = FitFunction->GetNDF();
383  for (int i = 0; i < 4; i++) {
384  par[i] = FitFunction->GetParameter(i);
385  errpar[i] = FitFunction->GetParError(i);
386  }
387 
388  } else { // The Selected Histo has NOT been fitted. Fit Now
389  double m = h->GetMean();
390  double r = h->GetRMS();
391  double maxval = h->GetMaximum();
392  std::array<double, 4> adcpar{};
393  adcpar[0] = maxval * 2.;
394  adcpar[1] = m;
395  adcpar[2] = r;
396  adcpar[3] = r / 3.;
397 
398  std::unique_ptr<TF1> AdcSpectrum = std::make_unique<TF1>(
399  "AdcSpectrum", " ([0]*exp((x-[1])/[2]))/ (1.+exp((x-[1])/[3])) ", m_settings->minAdc(), m_settings->maxAdc());
400  AdcSpectrum->SetParameters(adcpar.data());
401  double fitMin = m - (3 * r);
402  double fitMax = m + (3 * r);
403  if (fitMin < adcThreshold) fitMin = adcThreshold;
404  if (fitMax > 300) fitMax = 300.;
405  h->Fit("AdcSpectrum", "Q", " ", fitMin, fitMax);
406  chi2 = AdcSpectrum->GetChisquare();
407  ndof = AdcSpectrum->GetNDF();
408  for (int i = 0; i < 4; i++) {
409  par[i] = AdcSpectrum->GetParameter(i);
410  errpar[i] = AdcSpectrum->GetParError(i);
411  }
412  if (log.level() <= MSG::VERBOSE)
413  log << MSG::VERBOSE << "chi2/ndof=" << chi2 / ndof << " "
414  << "Mean=" << m << " "
415  << "RMS=" << r << " par 0 1 2 3 " << par[0] << " " << par[1] << " " << par[2] << " " << par[3] << endmsg;
416  }
417  // THE NEW HISTOGRAM HAS BEEN FITTED
418  }
419 
420  stc.adcCal = par[1];
421  fi.adc_par[0] = par[1];
422  fi.adc_par[1] = par[2];
423  fi.adc_err[0] = errpar[1];
424  fi.adc_err[1] = errpar[2];
425  if (std::abs(ndof) == 0) {
426  fi.adc_chi2 = -1;
427  } else {
428  fi.adc_chi2 = chi2 / static_cast<float>(ndof);
429  }
430  }
431  }
432 
433  void T0CalibrationClassic::searchParams(TH1 *h, double *p, int np) {
434  int nbinsX = h->GetNbinsX();
435  double sizeX = h->GetBinWidth(1);
436  double oldSizeX = sizeX;
437  int RebinFactor = static_cast<int>(10. / sizeX);
438  // extract starting values for fit params p[np] from the Time Spectrum h
439  std::unique_ptr<TH1> hnew{h->Rebin(RebinFactor, "hnew")}; // creates a new histogram hnew
440  // merging 5 bins of h1 in one bin
441  MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
442  if (log.level() <= MSG::DEBUG)
443  log << MSG::DEBUG << "nbinsx,sizex,rebinfactor=" << nbinsX << " " << sizeX << " " << RebinFactor << endmsg;
444  float minDeriv(9999.);
445  int minDerivBin(0);
446  sizeX = hnew->GetBinWidth(1);
447  int newbins = hnew->GetNbinsX();
448  for (int i = 0; i < np; ++i) {
449  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "i,p(i) " << i << " " << p[i] << endmsg;
450  }
451  for (int i = 0; i < newbins - 1; ++i) {
452  if (hnew->GetBinContent(i) - hnew->GetBinContent(i + 1) < minDeriv) {
453  minDeriv = hnew->GetBinContent(i) - hnew->GetBinContent(i + 1);
454  minDerivBin = i;
455  }
456  }
457  float t0guess = hnew->GetBinCenter(minDerivBin);
458  if (minDerivBin < newbins - 1) { t0guess += (hnew->GetBinCenter(minDerivBin + 1) - hnew->GetBinCenter(minDerivBin)) / 2.; }
459  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " t0guess is " << t0guess << endmsg;
460  //
461  // =================== Noise level search ===================================
462  //
463  float noise(0);
464  int numOfBins(10), numOfBinsOffset(3);
465  int imin, imax;
466  if (minDerivBin > numOfBins + numOfBinsOffset) {
467  imin = minDerivBin - numOfBins - numOfBinsOffset;
468  imax = minDerivBin - numOfBinsOffset;
469  } else {
470  imin = 0;
471  if (minDerivBin > numOfBinsOffset) {
472  imax = minDerivBin - numOfBinsOffset;
473  } else {
474  imax = minDerivBin;
475  }
476  }
477  int icount(0);
478  for (int i = imin; i <= imax; ++i) {
479  noise += hnew->GetBinContent(i);
480  icount++;
481  }
482 
483  noise = noise / (float)(icount);
484  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " noise is " << noise << endmsg;
485  //
486  // =================== Normalization =========================================
487  //
488  int t0bin = minDerivBin;
489  int ix1 = t0bin + (int)(50 / sizeX);
490  int ix2 = t0bin + (int)(500 / sizeX);
491  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "t0bin,ix1,ix2 " << t0bin << " " << ix1 << " " << ix2 << endmsg;
492  float P1 = p[1];
493  float P2 = p[2];
494  float P3 = p[3];
495  if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "P1,P2,P3 start are " << P1 << " " << P2 << " " << P3 << endmsg;
496  p[0] = noise;
497  p[4] = t0guess;
498  p[5] = p[4] + 700;
499  p[1] = 20.;
500  p[2] = 10.;
501  if (0 < ix1 && ix1 < newbins && 0 < ix2 && ix2 < newbins) {
502  float a1 = hnew->GetBinCenter(ix1);
503  float a2 = hnew->GetBinCenter(ix2);
504  float cont1 = hnew->GetBinContent(ix1);
505  float cont2 = hnew->GetBinContent(ix2);
506  if (cont1 > 0. && cont2 > 0.) {
507  float A1 = std::exp(-(a1 - t0guess) / P3);
508  float A2 = std::exp(-(a2 - t0guess) / P3);
509  // do not forget rebinning!
510  P2 = (cont1 / cont2 - 1.) / (A1 - cont1 / cont2 * A2);
511  P1 = cont1 / (1 + P2 * A1);
512  P1 = P1 * oldSizeX / sizeX;
513  P2 = P2 * oldSizeX / sizeX;
514  if (log.level() <= MSG::DEBUG) {
515  log << MSG::DEBUG << "a1,a2 " << a1 << " " << a2 << " cont1,cont2 " << cont1 << " " << cont2 << " A1,A2 " << A1 << " "
516  << A2 << endmsg;
517  log << MSG::DEBUG << " t0Guess .... P1, P2 " << P1 << " " << P2 << endmsg;
518  }
519  p[1] = P1;
520  p[2] = P2;
521  }
522  }
523  return;
524  }
525 
526  IMdtCalibration::MdtCalibOutputPtr T0CalibrationClassic::getResults() const {
527  return std::make_shared<T0CalibrationOutput>(m_result.get());
528  }
529 
530  T0ClassicHistos *T0CalibrationClassic::getHistos(unsigned int idtube) {
531  std::string HistoId;
532  if ((int)(idtube / 100000000) == 9) {
533  int mezz = (idtube) % (900000000);
534  HistoId = std::format("time_mezz_{:}", mezz);
535  } else if (idtube == 0) {
536  HistoId = "time";
537  } else if (idtube == 1) {
538  HistoId = "time_ML1";
539  } else if (idtube == 2) {
540  HistoId = "time_ML2";
541  } else if (idtube == 11) {
542  HistoId = "time_ML1_series1";
543  } else if (idtube == 12) {
544  HistoId = "time_ML1_series2";
545  } else if (idtube == 13) {
546  HistoId = "time_ML1_series3";
547  } else if (idtube == 21) {
548  HistoId = "time_ML2_series1";
549  } else if (idtube == 22) {
550  HistoId = "time_ML2_series2";
551  } else if (idtube == 23) {
552  HistoId = "time_ML2_series3";
553  } else {
554  MuonFixedId fId(idtube);
555  int nML = fId.mdtMultilayer();
556  int nL = fId.mdtTubeLayer();
557  int nT = fId.mdtTube();
558  int tubeid = (nML * 1000) + (nL * 100) + nT;
559  int nstat = fId.stationName();
560  std::string stationName = fId.stationNumberToFixedStationString(nstat);
561  int eta = fId.eta();
562  int phi = fId.phi();
563  HistoId = std::format("time_{:}_{:}_{:}_{:}", stationName, eta, phi, tubeid);
564  }
565  T0ClassicHistos *ret = nullptr;
566  TH1F *timeHis = (TH1F *)m_regiondir->Get(HistoId.c_str());
567  // We will either find the pointer to an existing histogram or create a new one
568  if (!timeHis) { // book histo if it does not exist
569  ret = bookHistos(idtube);
570  } else { // else loop over m_histos histograms to look for the set of histos of tube idtube
571  for (const std::unique_ptr<T0ClassicHistos> &hist : m_histos) {
572  if (hist->time.get() == timeHis) {
573  ret = hist.get();
574  break;
575  }
576  }
577  }
578  return ret;
579  } // end T0CalibrationClassic::getHistos
580 
581  T0ClassicHistos *T0CalibrationClassic::bookHistos(unsigned int idtube) {
582  std::unique_ptr<T0ClassicHistos> histos = std::make_unique<T0ClassicHistos>();
583  std::string histonametdc{}, histonameadc{};
584 
585  histos->id = idtube;
586  m_regiondir->cd();
587  if ((int)(idtube / 100000000) == 9) {
588  int mezz = (idtube) % (900000000);
589  histonametdc = std::format("time_mezz_{:}", mezz);
590  histonameadc = std::format("charge_mezz_{:}", mezz);
591  } else if (idtube == 0) {
592  histonametdc = "time";
593  histonameadc = "charge";
594  } else if (idtube == 1) {
595  histonametdc = "time_ML1";
596  histonameadc = "charge_ML1";
597  } else if (idtube == 2) {
598  histonametdc = "time_ML2";
599  histonameadc = "charge_ML2";
600  } else if (idtube == 11) {
601  histonametdc = "time_ML1_series1";
602  histonameadc = "charge_ML1_series1";
603  } else if (idtube == 12) {
604  histonametdc = "time_ML1_series2";
605  histonameadc = "charge_ML1_series2";
606  } else if (idtube == 13) {
607  histonametdc = "time_ML1_series3";
608  histonameadc = "charge_ML1_series3";
609  } else if (idtube == 21) {
610  histonametdc = "time_ML2_series1";
611  histonameadc = "charge_ML2_series1";
612  } else if (idtube == 22) {
613  histonametdc = "time_ML2_series2";
614  histonameadc = "charge_ML2_series2";
615  } else if (idtube == 23) {
616  histonametdc = "time_ML2_series3";
617  histonameadc = "charge_ML2_series3";
618  } else {
619  MuonFixedId fId(idtube);
620  int nML = fId.mdtMultilayer();
621  int nL = fId.mdtTubeLayer();
622  int nT = fId.mdtTube();
623  int nstat = fId.stationName();
624  int tubeid = (nML * 1000) + (nL * 100) + nT;
625  std::string stationName = fId.stationNumberToFixedStationString(nstat);
626  int eta = fId.eta();
627  int phi = fId.phi();
628  histonametdc = std::format("time_{:}_{:}_{:}_{:}", stationName, eta, phi, tubeid);
629  histonameadc = std::format("charge_{:}_{:}_{:}_{:}", stationName, eta, phi, tubeid);
630  }
631 
632  histos->time =
633  std::make_unique<TH1F>(histonametdc.c_str(), "Drift Time", m_settings->binTime(), m_settings->minTime(), m_settings->maxTime());
634  histos->adc = std::make_unique<TH1F>(histonameadc.c_str(), "ADC", m_settings->binAdc(), m_settings->minAdc(), m_settings->maxAdc());
635 
636  m_histos.emplace_back(std::move(histos));
637  return m_histos.back().get();
638  }
639 
640  void T0CalibrationClassic::setInput(const IMdtCalibrationOutput *calib_in) {
641  // This method is called both by the event loop and by the tool.
642  // Only the call from the tool is relevant for this implementation
643  // and should be performed only once.
644 
645  if (m_result) return;
646 
647  const T0CalibrationOutput *t0Input = dynamic_cast<const T0CalibrationOutput *>(calib_in);
648  if (t0Input) {
649  m_result = std::make_unique<MdtTubeFitContainer>(*t0Input->t0s());
650  m_result->setImplementation("T0CalibrationClassic");
651  }
652  }
653 
654  bool T0CalibrationClassic::converged() const { return m_converged; }
655 
656 } // namespace MuonCalib
MuonCalib::T0CalibrationClassic::bookHistos
T0ClassicHistos * bookHistos(unsigned int idtube)
booking of histograms
Definition: T0CalibrationClassic.cxx:581
MuonCalib::T0ClassicSettings::numParams
int numParams() const
Definition: T0CalibrationClassic.h:69
MuonCalib::T0CalibrationClassic::m_currentItnum
int m_currentItnum
current iteration (always 1?)
Definition: T0CalibrationClassic.h:148
MuonCalib::MdtTubeCalibContainer::SingleTubeCalib::t0
float t0
< relative t0 in chamber (ns)
Definition: MdtTubeCalibContainer.h:21
MuonCalib::MdtTubeCalibContainer::SingleTubeCalib::statusCode
uint8_t statusCode
Definition: MdtTubeCalibContainer.h:25
MuonCalib::MuonFixedId::isValid
bool isValid() const
check validity of the identifier.
Definition: MuonFixedId.h:555
beamspotman.r
def r
Definition: beamspotman.py:676
IDTPM::ndof
float ndof(const U &p)
Definition: TrackParametersHelper.h:134
MuonCalib::IMdtCalibration::MdtCalibOutputPtr
std::shared_ptr< IMdtCalibrationOutput > MdtCalibOutputPtr
Definition: IMdtCalibration.h:30
MuonCalib::T0CalibrationClassic::m_result
std::unique_ptr< MdtTubeFitContainer > m_result
tube constants
Definition: T0CalibrationClassic.h:153
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
MuonCalibSegment.h
vtune_athena.format
format
Definition: vtune_athena.py:14
MuonCalib::T0CalibrationClassic::m_file
std::unique_ptr< TFile > m_file
pointer to the histogram file
Definition: T0CalibrationClassic.h:149
MuonCalib::MdtTubeFitContainer::SingleTubeFit
Definition: MdtTubeFitContainer.h:18
MuonCalib::T0CalibrationClassic::m_name
std::string m_name
calibration region name
Definition: T0CalibrationClassic.h:147
dumpTgcDigiDeadChambers.stationName
dictionary stationName
Definition: dumpTgcDigiDeadChambers.py:30
MuonCalib::TimeSpectrum_func
Double_t TimeSpectrum_func(Double_t *xx, Double_t *par)
Definition: T0CalibrationClassic.cxx:24
plotmaker.hist
hist
Definition: plotmaker.py:148
MuonCalib::T0CalibrationClassic::setInput
void setInput(const IMdtCalibrationOutput *input) override
unused
Definition: T0CalibrationClassic.cxx:640
MuonCalib::T0CalibrationClassic::getResults
virtual MdtCalibOutputPtr getResults() const override
Definition: T0CalibrationClassic.cxx:526
MuonCalib::MuonFixedId::is_mdt
bool is_mdt() const
Definition: MuonFixedId.h:559
MuonCalib::IMdtCalibration
Definition: IMdtCalibration.h:25
MuonCalib::MuonFixedId::mdtMultilayer
int mdtMultilayer() const
Mdt specific:
Definition: MuonFixedId.h:835
MuonCalib::MuonCalibSegment
Definition: MuonCalibSegment.h:39
T0CalibrationClassic.h
MuonCalib::T0ClassicSettings::binAdc
int binAdc() const
Definition: T0CalibrationClassic.h:54
MuonCalib::T0CalibrationClassic::m_regiondir
TDirectory * m_regiondir
pointer to the ROOT directory
Definition: T0CalibrationClassic.h:150
PlotPulseshapeFromCool.np
np
Definition: PlotPulseshapeFromCool.py:64
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
x
#define x
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition: getMessageSvc.cxx:20
MuonCalib::T0ClassicSettings::maxAdc
double maxAdc() const
Definition: T0CalibrationClassic.h:52
MuonCalib::T0CalibrationClassic::m_delete_settings
bool m_delete_settings
Definition: T0CalibrationClassic.h:154
MuonCalib::T0ClassicSettings::printLevel
int printLevel() const
Definition: T0CalibrationClassic.h:75
MuonCalib::T0CalibrationOutput
Definition: T0CalibrationOutput.h:19
MuonCalib::MuonFixedId::stationNameString
std::string stationNameString() const
Definition: MuonFixedId.h:655
MuonCalib::T0ClassicHistos
Definition: T0CalibrationClassic.h:108
T0CalibrationOutput.h
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
MuonCalib::T0CalibrationClassic::handleSegment
bool handleSegment(MuonCalibSegment &seg)
fill tube spectra
Definition: T0CalibrationClassic.cxx:62
MuonCalib::T0ClassicHistos::time
std::unique_ptr< TH1 > time
time spectrum
Definition: T0CalibrationClassic.h:110
lumiFormat.i
int i
Definition: lumiFormat.py:85
xmin
double xmin
Definition: listroot.cxx:60
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
python.LArMinBiasAlgConfig.int
int
Definition: LArMinBiasAlgConfig.py:59
MuonCalib::T0CalibrationClassic::T0CalibrationClassic
T0CalibrationClassic(const std::string &name, const T0ClassicSettings *settings)
constructor
Definition: T0CalibrationClassic.cxx:30
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:525
MuonCalib::T0ClassicSettings::entries
int entries() const
Definition: T0CalibrationClassic.h:64
MdtIdHelper
Definition: MdtIdHelper.h:61
MuonCalib
CscCalcPed - algorithm that finds the Cathode Strip Chamber pedestals from an RDO.
Definition: CscCalcPed.cxx:22
MuonCalib::T0ClassicSettings::print
void print() const
a method to dump all the settings
Definition: T0CalibrationClassic.h:77
MuonCalib::T0ClassicSettings::minTime
double minTime() const
Definition: T0CalibrationClassic.h:56
MuonCalib::T0CalibrationClassic::m_settings
const T0ClassicSettings * m_settings
pointer to the settings
Definition: T0CalibrationClassic.h:145
MuonCalib::T0CalibrationClassic::converged
bool converged() const
return m_converged (always false?)
Definition: T0CalibrationClassic.cxx:654
MuonCalib::MuonFixedId::stationNumberToFixedStationString
static std::string stationNumberToFixedStationString(const int station)
Definition: MuonFixedId.h:747
imax
int imax(int i, int j)
Definition: TileLaserTimingTool.cxx:33
MuonCalib::T0CalibrationClassic::m_converged
bool m_converged
convergence status
Definition: T0CalibrationClassic.h:146
MuonCalib::IMdtCalibration::name
virtual std::string name() const
returns name (region) of instance
Definition: IMdtCalibration.h:49
MuonCalib::MdtTubeCalibContainer::SingleTubeCalib::adcCal
float adcCal
quality flag for the SingleTubeCalib constants: 0 all ok, 1 no hits found, 2 too few hits,...
Definition: MdtTubeCalibContainer.h:23
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
MuonCalib::MuonFixedId::stationName
int stationName() const
Definition: MuonFixedId.h:651
MuonCalib::T0CalibrationOutput::t0s
const MdtTubeFitContainer * t0s() const
Definition: T0CalibrationOutput.h:27
MuonCalib::T0ClassicSettings::params
const std::array< double, 8 > & params() const
Definition: T0CalibrationClassic.h:71
MuonCalib::MuonFixedId
Definition: MuonFixedId.h:50
MuonCalib::IMdtCalibrationOutput
Definition: IMdtCalibrationOutput.h:28
find_data.full
full
Definition: find_data.py:27
MuonCalib::MuonFixedId::phi
int phi() const
Definition: MuonFixedId.h:704
MuonCalib::T0ClassicSettings::maxTime
double maxTime() const
Definition: T0CalibrationClassic.h:58
python.testIfMatch.matrix
matrix
Definition: testIfMatch.py:66
h
DEBUG
#define DEBUG
Definition: page_access.h:11
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
xmax
double xmax
Definition: listroot.cxx:61
MuonCalib::T0ClassicSettings
Definition: T0CalibrationClassic.h:30
entries
double entries
Definition: listroot.cxx:49
MuonCalib::T0CalibrationClassic::getHistos
T0ClassicHistos * getHistos(unsigned int idtube)
retrieve pointer for tube idtube histograms
Definition: T0CalibrationClassic.cxx:530
MuonCalib::T0CalibrationClassic::doTimeFit
void doTimeFit(T0ClassicHistos &, MdtTubeFitContainer::SingleTubeFit &, MdtTubeFitContainer::SingleTubeCalib &)
fit time spectrum
Definition: T0CalibrationClassic.cxx:159
checkFileSG.fi
fi
Definition: checkFileSG.py:65
MuonCalib::MuonFixedId::mdtMezzanine
int mdtMezzanine() const
Mdt specific: compute the mezzanine number.
Definition: MuonFixedId.h:779
MuonCalib::MuonFixedId::mdtTube
int mdtTube() const
Mdt specific:
Definition: MuonFixedId.h:775
MuonCalib::MuonFixedId::eta
int eta() const
Definition: MuonFixedId.h:681
MuonCalib::MdtTubeCalibContainer::SingleTubeCalib
Definition: MdtTubeCalibContainer.h:19
MuonCalib::IMdtCalibration::MuonSegVec
std::vector< std::shared_ptr< MuonCalibSegment > > MuonSegVec
Definition: IMdtCalibration.h:27
checkCorrelInHIST.histos
dictionary histos
Definition: checkCorrelInHIST.py:413
MuonCalib::T0ClassicSettings::binTime
int binTime() const
Definition: T0CalibrationClassic.h:60
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
MuonCalib::MuonCalibSegment::mdtHOT
const MdtHitVec & mdtHOT() const
retrieve the full set of MdtCalibHitBase s assigned to this segment
Definition: MuonCalibSegment.cxx:147
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:13
PowhegControl_ttFCNC_NLO.params
params
Definition: PowhegControl_ttFCNC_NLO.py:226
MuonCalib::T0ClassicSettings::fitTime
bool fitTime() const
Definition: T0CalibrationClassic.h:62
MuonCalib::T0CalibrationClassic::analyse
bool analyse()
extract parameters from spectra
Definition: T0CalibrationClassic.cxx:119
MuonCalib::T0ClassicHistos::adc
std::unique_ptr< TH1 > adc
adc spectrum
Definition: T0CalibrationClassic.h:111
MuonCalib::T0CalibrationClassic::m_histos
std::vector< std::unique_ptr< T0ClassicHistos > > m_histos
vector of pointers to tube histograms
Definition: T0CalibrationClassic.h:152
MuonCalib::T0CalibrationClassic::~T0CalibrationClassic
~T0CalibrationClassic()
destructor
Definition: T0CalibrationClassic.cxx:56
MuonCalib::T0CalibrationClassic::searchParams
void searchParams(TH1 *h, double *p, int np)
estimate initial parameters for time spectrum fit from the spectrum itself
Definition: T0CalibrationClassic.cxx:433
MuonCalib::T0CalibrationClassic::analyseSegments
virtual MdtCalibOutputPtr analyseSegments(const MuonSegVec &segs) override
new interface function
Definition: T0CalibrationClassic.cxx:153
MuonFixedId.h
MuonCalib::T0ClassicSettings::minAdc
double minAdc() const
Access functions to values of private settings.
Definition: T0CalibrationClassic.h:50
WriteCellNoiseToCool.noise
noise
Definition: WriteCellNoiseToCool.py:380
MuonCalib::MuonFixedId::mdtTubeLayer
int mdtTubeLayer() const
Mdt specific:
Definition: MuonFixedId.h:813
MuonCalib::MuonCalibSegment::MdtHitPtr
std::shared_ptr< MdtCalibHitBase > MdtHitPtr
typedef for a collection of MdtCalibHitBase s
Definition: MuonCalibSegment.h:44
python.LArMinBiasAlgConfig.float
float
Definition: LArMinBiasAlgConfig.py:65
MuonCalib::T0ClassicHistos::id
int id
tube identifier
Definition: T0CalibrationClassic.h:113
MuonCalib::T0CalibrationClassic::doAdcFit
void doAdcFit(T0ClassicHistos &, MdtTubeFitContainer::SingleTubeFit &, MdtTubeFitContainer::SingleTubeCalib &)
fit adc spectrum
Definition: T0CalibrationClassic.cxx:327
Identifier
Definition: IdentifierFieldParser.cxx:14
Generated on Sun Mar 23 2025 21:19:13 for ATLAS Offline Software by   doxygen 1.8.18