ATLAS Offline Software
Loading...
Searching...
No Matches
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
20namespace 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(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
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
endmsg
#define endmsg
Definition AnalysisConfig_Ntuple.cxx:54
hit
bool hit(const Container &ids, int pdgId)
Definition JetIRCSafeLabelTool.cxx:64
imax
int imax(int i, int j)
Definition TileLaserTimingTool.cxx:33
x
#define x
h
Header file for AthHistogramAlgorithm.
MdtIdHelper::channelID
Identifier channelID(int stationName, int stationEta, int stationPhi, int multilayer, int tubeLayer, int tube) const
Definition MdtIdHelper.cxx:663
MuonCalib::IMdtCalibrationOutput
Interface to pass calibration output during calibration.
Definition IMdtCalibrationOutput.h:28
MuonCalib::IMdtCalibration::IMdtCalibration
IMdtCalibration(const std::string &name)
constructor, string used to identify the instance
Definition IMdtCalibration.h:34
MuonCalib::IMdtCalibration::name
virtual std::string name() const
returns name (region) of instance
Definition IMdtCalibration.h:49
MuonCalib::IMdtCalibration::MdtCalibOutputPtr
std::shared_ptr< IMdtCalibrationOutput > MdtCalibOutputPtr
Definition IMdtCalibration.h:30
MuonCalib::IMdtCalibration::MuonSegVec
std::vector< std::shared_ptr< MuonCalibSegment > > MuonSegVec
Definition IMdtCalibration.h:27
MuonCalib::MuonCalibSegment
A MuonCalibSegment is a reconstructed three dimensional track segment in the MuonSpectrometer.
Definition MuonCalibSegment.h:39
MuonCalib::MuonCalibSegment::MdtHitPtr
std::shared_ptr< MdtCalibHitBase > MdtHitPtr
typedef for a collection of MdtCalibHitBase s
Definition MuonCalibSegment.h:43
MuonCalib::MuonCalibSegment::mdtHOT
const MdtHitVec & mdtHOT() const
retrieve the full set of MdtCalibHitBase s assigned to this segment
Definition MuonCalibSegment.cxx:147
MuonCalib::MuonFixedId
Implements fixed identifiers not dependent upon Athena Identifier for internal use in the calibration...
Definition MuonFixedId.h:50
MuonCalib::MuonFixedId::mdtMezzanine
int mdtMezzanine() const
Mdt specific: compute the mezzanine number.
Definition MuonFixedId.h:779
MuonCalib::MuonFixedId::stationNumberToFixedStationString
static std::string stationNumberToFixedStationString(const int station)
Definition MuonFixedId.h:747
MuonCalib::MuonFixedId::mdtTubeLayer
int mdtTubeLayer() const
Mdt specific:
Definition MuonFixedId.h:814
MuonCalib::MuonFixedId::mdtTube
int mdtTube() const
Mdt specific:
Definition MuonFixedId.h:775
MuonCalib::MuonFixedId::isValid
bool isValid() const
check validity of the identifier.
Definition MuonFixedId.h:555
MuonCalib::MuonFixedId::stationNameString
std::string stationNameString() const
Definition MuonFixedId.h:655
MuonCalib::MuonFixedId::mdtMultilayer
int mdtMultilayer() const
Mdt specific:
Definition MuonFixedId.h:836
MuonCalib::T0CalibrationClassic::m_name
std::string m_name
calibration region name
Definition T0CalibrationClassic.h:147
MuonCalib::T0CalibrationClassic::m_currentItnum
int m_currentItnum
current iteration (always 1?)
Definition T0CalibrationClassic.h:148
MuonCalib::T0CalibrationClassic::m_file
std::unique_ptr< TFile > m_file
pointer to the histogram file
Definition T0CalibrationClassic.h:149
MuonCalib::T0CalibrationClassic::analyseSegments
virtual MdtCalibOutputPtr analyseSegments(const MuonSegVec &segs) override
new interface function
Definition T0CalibrationClassic.cxx:153
MuonCalib::T0CalibrationClassic::handleSegment
bool handleSegment(MuonCalibSegment &seg)
fill tube spectra
Definition T0CalibrationClassic.cxx:62
MuonCalib::T0CalibrationClassic::getResults
virtual MdtCalibOutputPtr getResults() const override
Definition T0CalibrationClassic.cxx:526
MuonCalib::T0CalibrationClassic::setInput
void setInput(const IMdtCalibrationOutput *input) override
unused
Definition T0CalibrationClassic.cxx:640
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::m_settings
const T0ClassicSettings * m_settings
pointer to the settings
Definition T0CalibrationClassic.h:145
MuonCalib::T0CalibrationClassic::m_result
std::unique_ptr< MdtTubeFitContainer > m_result
tube constants
Definition T0CalibrationClassic.h:153
MuonCalib::T0CalibrationClassic::converged
bool converged() const
return m_converged (always false?)
Definition T0CalibrationClassic.cxx:654
MuonCalib::T0CalibrationClassic::m_regiondir
TDirectory * m_regiondir
pointer to the ROOT directory
Definition T0CalibrationClassic.h:150
MuonCalib::T0CalibrationClassic::m_histos
std::vector< std::unique_ptr< T0ClassicHistos > > m_histos
vector of pointers to tube histograms
Definition T0CalibrationClassic.h:152
MuonCalib::T0CalibrationClassic::doAdcFit
void doAdcFit(T0ClassicHistos &, MdtTubeFitContainer::SingleTubeFit &, MdtTubeFitContainer::SingleTubeCalib &)
fit adc spectrum
Definition T0CalibrationClassic.cxx:327
MuonCalib::T0CalibrationClassic::doTimeFit
void doTimeFit(T0ClassicHistos &, MdtTubeFitContainer::SingleTubeFit &, MdtTubeFitContainer::SingleTubeCalib &)
fit time spectrum
Definition T0CalibrationClassic.cxx:159
MuonCalib::T0CalibrationClassic::getHistos
T0ClassicHistos * getHistos(unsigned int idtube)
retrieve pointer for tube idtube histograms
Definition T0CalibrationClassic.cxx:530
MuonCalib::T0CalibrationClassic::T0CalibrationClassic
T0CalibrationClassic(const std::string &name, const T0ClassicSettings *settings)
constructor
Definition T0CalibrationClassic.cxx:30
MuonCalib::T0CalibrationClassic::bookHistos
T0ClassicHistos * bookHistos(unsigned int idtube)
booking of histograms
Definition T0CalibrationClassic.cxx:581
MuonCalib::T0CalibrationClassic::analyse
bool analyse()
extract parameters from spectra
Definition T0CalibrationClassic.cxx:119
MuonCalib::T0CalibrationOutput
class for the communication of the results of T0 calibration algorithms
Definition T0CalibrationOutput.h:19
MuonCalib::T0CalibrationOutput::t0s
const MdtTubeFitContainer * t0s() const
Definition T0CalibrationOutput.h:27
MuonCalib::T0ClassicHistos
Tube histograms used in T0 calibration.
Definition T0CalibrationClassic.h:108
MuonCalib::T0ClassicHistos::adc
std::unique_ptr< TH1 > adc
adc spectrum
Definition T0CalibrationClassic.h:111
MuonCalib::T0ClassicHistos::time
std::unique_ptr< TH1 > time
time spectrum
Definition T0CalibrationClassic.h:110
MuonCalib::T0ClassicSettings
Settings for the T0 calibration (histogram booking and fitting parameters).
Definition T0CalibrationClassic.h:30
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition comparitor.cxx:525
r
int r
Definition globals.cxx:22
xmax
double xmax
Definition listroot.cxx:61
entries
double entries
Definition listroot.cxx:49
xmin
double xmin
Definition listroot.cxx:60
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition getMessageSvc.cxx:20
MuonCalib
CscCalcPed - algorithm that finds the Cathode Strip Chamber pedestals from an RDO.
Definition CscCalcPed.cxx:21
MuonCalib::TimeSpectrum_func
Double_t TimeSpectrum_func(Double_t *xx, Double_t *par)
Definition T0CalibrationClassic.cxx:24
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:28
MuonCalib::MdtTubeFitContainer::SingleTubeFit::statistics
int statistics
< number of hits used for the fit
Definition MdtTubeFitContainer.h:20
Generated on for ATLAS Offline Software by doxygen 1.17.0