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