ATLAS Offline Software
CscFit.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // a class to fit CSC strips in a plane into clusters //
7 // //
8 // BNL March 26 2003 --- Ketevi A. Assamagan //
10 
11 #include "CscFit.h"
12 
13 #include <cmath>
14 
15 #include "AthenaKernel/getMessageSvc.h"
16 #include "GaudiKernel/MsgStream.h"
17 #include "TF1.h"
18 #include "TH1F.h"
19 
20 CscFit::CscFit(double sigma) { m_sigma = sigma; }
21 
29 
30 void CscFit::cscfit(double* qstr, int& maxStrip, double& thr, double& da, int& ncl, double* sig, double* zpos, double& noise) const {
31 #ifndef NDEBUG
32  MsgStream log(Athena::getMessageSvc(), "CscFit");
33 #endif
34 
35  const int NN = 100;
36  int NPT;
37  int ii;
38  int imax;
39  int iquest1, iquest2;
40  int ithrL, ithrR, is; // i0=0;
41 
42  bool first = true;
43 
44  double zl, zr;
45  double q1, q2;
46  double /* ddz=0, */ zfcsc1 = -1000.0, zcscl = 0, zcscr = 0; // dlr=-1000.0;
47 #ifndef NDEBUG
48  double fErr1 = 10000.0;
49 #endif
50  double dzcut = da / 25.0;
51 
52  int L = maxStrip;
53 
54  ithrL = 1;
55  ithrR = 1;
56 
57  if (first) {
58  // i0 = L/2;
59  first = false;
60  }
61  q2 = 0;
62  ncl = 0;
63  for (is = 1; is < L; is++) {
64  q1 = q2;
65  q2 = *(qstr + is);
66  if (is == (L - 1)) q2 = 0;
67  if (q1 < thr && q2 >= thr) ithrL = is;
68  if (q1 < thr || q2 >= thr) continue;
69  ithrR = is - 1;
70 
72  imax = icmax(qstr, ithrL, std::min(L - 1, (ithrL + 2)));
73  iquest1 = std::max(ithrL - 2, 0);
74  iquest2 = imax;
75  double al = qstr[imax] / 2.0;
76  NPT = 0;
77  double sum = 0;
78  double norm = 0;
79  for (ii = 0; ii < L; ii++)
80  if (ii >= iquest1 && ii <= iquest2 && *(qstr + ii) > 1) NPT++;
81  double lowerBin = da * iquest1 - da * L / 2;
82  double higherBin = da * iquest2 - da * L / 2;
83  TH1F* hist = new TH1F("csc", "csc", NPT, lowerBin, higherBin);
84  for (ii = 0; ii < L; ii++) {
85  if (ii >= iquest1 && ii <= iquest2 && *(qstr + ii) > 1) {
86  double y = *(qstr + ii);
87  double x = da * ((ii + 0.5) - L / 2);
88  sum += y * x;
89  norm += y;
90  hist->Fill(x, y);
91  }
92  }
93  zl = sum / norm;
94 
95  TF1* func = new TF1("f1gauss", f1gauss, lowerBin, higherBin, 4);
96  if (NPT >= 4) {
97  Double_t par[4];
98  Double_t x0 = da * ((imax + 0.5) - L / 2);
99  Double_t s0 = da * m_sigma;
100  func->SetParameters(qstr[imax], x0, s0, 0.0);
101  int status = hist->Fit("f1gauss", "QRI");
102  func->GetParameters(&par[0]);
103 #ifndef NDEBUG
104  fErr1 = func->GetParError(1);
105  log << MSG::DEBUG << "One Gaussing Left Fit: "
106  << "Avg = " << zl << " Mean = " << par[1] << " Error = " << fErr1 << " status = " << status << endmsg;
107 #endif
108  if (status == 0) {
109  al = (double)par[0];
110  zl = (double)par[1];
111  }
112  }
113  delete func;
114  delete hist;
115 
117  imax = icmax(qstr, std::max(0, (ithrR - 2)), ithrR);
118  iquest1 = imax;
119  iquest2 = std::min(ithrR + 2, L - 1);
120  double ar = qstr[imax] / 2.0;
121  NPT = 0;
122  sum = 0;
123  norm = 0;
124  for (ii = 0; ii < L; ii++)
125  if (ii >= iquest1 && ii <= iquest2 && *(qstr + ii) > 1) NPT++;
126  lowerBin = da * iquest1 - da * L / 2;
127  higherBin = da * iquest2 - da * L / 2;
128  hist = new TH1F("csc", "csc", NPT, lowerBin, higherBin);
129  for (ii = 0; ii < L; ii++) {
130  if (ii >= iquest1 && ii <= iquest2 && *(qstr + ii) > 1) {
131  double y = *(qstr + ii);
132  double x = da * ((ii + 0.5) - L / 2);
133  sum += y * x;
134  norm += y;
135  hist->Fill(x, y);
136  }
137  }
138 
139  zr = sum / norm;
140 
141  func = new TF1("f1gauss", f1gauss, lowerBin, higherBin, 4);
142  if (NPT >= 4) {
143  Double_t par[4];
144  Double_t x0 = da * ((imax + 0.5) - L / 2);
145  Double_t s0 = da * m_sigma;
146  func->SetParameters(qstr[imax], x0, s0, 0.0);
147  int status = hist->Fit("f1gauss", "QRI");
148  func->GetParameters(&par[0]);
149 #ifndef NDEBUG
150  fErr1 = func->GetParError(1);
151  log << MSG::DEBUG << "One Gaussing Right Fit: "
152  << "Avg = " << zr << " Mean = " << par[1] << " Error = " << fErr1 << " status = " << status << endmsg;
153 #endif
154  if (status == 0) {
155  ar = (double)par[0];
156  zr = (double)par[1];
157  }
158  }
159  delete func;
160  delete hist;
161 
162  // ddz = fabs(zr-zl);
163 
165  iquest1 = std::max(ithrL - 2, 0);
166  iquest2 = std::min(ithrR + 2, L - 1);
167  int jmax = icmax(qstr, iquest1, iquest2) + 1;
168 
169  NPT = 0;
170  sum = 0;
171  norm = 0;
172  for (ii = 0; ii < L; ii++)
173  if (ii >= iquest1 && ii <= iquest2 && *(qstr + ii) > 1) NPT++;
174  lowerBin = da * iquest1 - da * L / 2;
175  higherBin = da * iquest2 - da * L / 2;
176  hist = new TH1F("csc", "csc", NPT, lowerBin, higherBin);
177  for (ii = 0; ii < L; ii++) {
178  if (ii >= iquest1 && ii <= iquest2 && *(qstr + ii) > 1) {
179  double y = *(qstr + ii);
180  double x = da * ((ii + 0.5) - L / 2);
181  hist->Fill(x, y);
182  sum += y * x;
183  norm += y;
184  }
185  }
186 
187  double z0 = sum / norm;
188  zfcsc1 = z0;
189  double sigMainFit = 0.0;
190 
191  bool mainFit = false;
192  func = new TF1("f1gauss", f1gauss, lowerBin, higherBin, 4);
193  if (NPT >= 4) {
194  Double_t par[4];
195  Double_t s0 = da * m_sigma;
196  sigMainFit = s0;
197  func->SetParameters(qstr[jmax], z0, s0, 0.0);
198  // func->FixParameter(3,s0);
199  int status = hist->Fit("f1gauss", "QRBI");
200  func->GetParameters(&par[0]);
201  if (status == 0) {
202  zfcsc1 = (double)par[1];
203  sigMainFit = (double)par[2];
204 #ifndef NDEBUG
205  fErr1 = func->GetParError(1);
206 #endif
207  mainFit = true;
208  }
209 #ifndef NDEBUG
210  log << MSG::DEBUG << "One Gaussing Main Fit: "
211  << "Avg = " << z0 << " Mean = " << par[1] << " Error = " << fErr1 << " baseline = " << par[3] << " status = " << status
212  << endmsg;
213 #endif
214  }
215  delete func;
216  delete hist;
217 
219  iquest1 = std::max(ithrL - 2, 0);
220  iquest2 = std::min(ithrR + 2, L - 1);
221  NPT = 0;
222  sum = 0;
223  norm = 0;
224  int max_strip = -1;
225  double max_charge = -1;
226  for (ii = 0; ii < L; ii++)
227  if (ii >= iquest1 && ii <= iquest2 && *(qstr + ii) > 1) NPT++;
228  lowerBin = da * iquest1 - da * L / 2;
229  higherBin = da * iquest2 - da * L / 2;
230  hist = new TH1F("csc", "csc", NPT, lowerBin, higherBin);
231  for (ii = 0; ii < L; ii++) {
232  if (ii >= iquest1 && ii <= iquest2 && *(qstr + ii) > 1) {
233  double y = *(qstr + ii);
234  double x = da * ((ii + 0.5) - L / 2);
235  sum += y * x;
236  norm += y;
237  if (y > max_charge) {
238  max_charge = y;
239  max_strip = ii;
240  }
241  hist->Fill(x, y);
242  }
243  }
244 
245  bool doubleGaussianFit = false;
246  func = new TF1("f2gauss", f2gauss, lowerBin, higherBin, 7);
247  if (NPT >= 7) {
248  Double_t par[7];
249  double_t s0 = da * m_sigma / 2.0;
250  if (mainFit) s0 = sigMainFit / 2.0;
251  func->SetParameters(al, zl, s0, ar, zr, s0, 0.0);
252  int status = hist->Fit("f2gauss", "QRI");
253  func->GetParameters(&par[0]);
254  if (status == 0) {
255  doubleGaussianFit = true;
256  zcscl = (double)par[1];
257  zcscr = (double)par[4];
258  }
259 #ifndef NDEBUG
260  log << MSG::DEBUG << "Double Gaussing Fit: "
261  << "Avg = " << z0 << " Mean One Gaussian = " << zfcsc1 << " Mean1 = " << par[1] << " Error1 = " << func->GetParError(1)
262  << " Mean2 = " << par[4] << " Error2 = " << func->GetParError(4) << " baseline = " << par[6] << " status = " << status
263  << endmsg;
264 #endif
265  }
266  delete func;
267  delete hist;
268 
269  double ysum = 0;
270  double ysum2 = 0;
271  int count = 0;
272  for (int i = max_strip - 1; i <= max_strip + 1; i++) {
273  if (i >= 0 && i < L) {
274  ysum += *(qstr + i);
275  count++;
276  ysum2 += count * count;
277  }
278  }
279 
281  double wErr = da / sqrt(12.0);
282  if (count > 1 && ysum > 0) wErr = noise * (da / (ysum)) * sqrt(2.0 * ysum2);
283 
286  double xp0 = 0, xp1 = 0, xp2 = 0;
287 
288  xp0 = zfcsc1;
289  xp1 = std::min(zcscl, zcscr);
290  xp2 = std::max(zcscl, zcscr);
291 
292  double dx0 = fabs(xp0 - z0);
293  double dx1 = fabs(xp1 - z0);
294  double dx2 = fabs(xp2 - z0);
295 
296  bool checkIt = (dx0 < dzcut) || (dx1 < dzcut) || (dx2 < dzcut);
297 
298  if (mainFit && doubleGaussianFit && checkIt) {
299 #ifndef NDEBUG
300  log << MSG::DEBUG << "Double Gaussian Fit converged ..." << endmsg;
301 #endif
302  ncl = std::min(NN, ncl + 1);
303  if (dx0 < dx1 && dx0 < dx2) {
304  *(zpos + ncl - 1) = xp0;
305  } else if (dx1 < dx2 && dx1 < dx0) {
306  *(zpos + ncl - 1) = xp1;
307  } else if (dx2 < dx0 && dx2 < dx1) {
308  *(zpos + ncl - 1) = xp2;
309  } else {
310  *(zpos + ncl - 1) = xp0;
311  }
312  *(sig + ncl - 1) = wErr;
313  } else if (mainFit && (dx0 < dzcut)) {
314 #ifndef NDEBUG
315  log << MSG::DEBUG << "One Gaussian Fit converged ..." << endmsg;
316 #endif
317  ncl = std::min(NN, ncl + 1);
318  *(zpos + ncl - 1) = xp0;
319  *(sig + ncl - 1) = wErr;
320  } else {
321  if (ysum > 0) {
322  ncl = std::min(NN, ncl + 1);
323 #ifndef NDEBUG
324  log << MSG::DEBUG << "Fit was not meaningful: using weighted average" << endmsg;
325 #endif
326  *(zpos + ncl - 1) = z0;
327  *(sig + ncl - 1) = wErr;
328  }
329  }
330  }
331 }
332 
334 int CscFit::icmax(const double* qstr, const int& i1, const int& i2) {
335  int i, ic;
336  double a;
337 
338  a = *(qstr + i1);
339  ic = i1;
340  for (i = i1 + 1; i <= i2; i++) {
341  if (*(qstr + i) < a) continue;
342  a = *(qstr + i);
343  ic = i;
344  }
345  return (ic);
346 }
347 
349 Double_t f1gauss(const Double_t* x, const Double_t* par) {
350  Double_t arg = 0;
351  arg = (x[0] - par[1]) / par[2];
352  Double_t fitval = par[3] + par[0] * std::exp(-0.5 * arg * arg);
353  return fitval;
354 }
355 
357 Double_t f2gauss(const Double_t* x, const Double_t* par) {
358  Double_t arg1 = 0;
359  Double_t arg2 = 0;
360  if (par[2] != 0) arg1 = (x[0] - par[1]) / par[2];
361  if (par[5] != 0) arg2 = (x[0] - par[4]) / par[5];
362  Double_t fitval = par[0] * std::exp(-0.5 * arg1 * arg1) + par[3] * std::exp(-0.5 * arg2 * arg2) + par[6];
363  return fitval;
364 }
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
max
#define max(a, b)
Definition: cfImp.cxx:41
getMessageSvc.h
singleton-like access to IMessageSvc via open function and helper
CscFit::icmax
static int icmax(const double *qstr, const int &i1, const int &i2)
a method to find a the index of the strip with the highest charge
Definition: CscFit.cxx:334
plotmaker.hist
hist
Definition: plotmaker.py:148
f2gauss
Double_t f2gauss(const Double_t *x, const Double_t *par)
two-Gaussian fitting method
Definition: CscFit.cxx:357
MCP::ScaleSmearParam::s0
@ s0
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
CscFit::m_sigma
double m_sigma
Definition: CscFit.h:35
x
#define x
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition: getMessageSvc.cxx:20
XMLtoHeader.count
count
Definition: XMLtoHeader.py:85
f1gauss
Double_t f1gauss(const Double_t *x, const Double_t *par)
one Gaussing fitting method
Definition: CscFit.cxx:349
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
CscFit::f2gauss
friend Double_t f2gauss(const Double_t *x, const Double_t *par)
two-Gaussian fitting method
Definition: CscFit.cxx:357
lumiFormat.i
int i
Definition: lumiFormat.py:92
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
python.BuildSignatureFlags.sig
sig
Definition: BuildSignatureFlags.py:215
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
CscFit::f1gauss
friend Double_t f1gauss(const Double_t *x, const Double_t *par)
one Gaussing fitting method
Definition: CscFit.cxx:349
imax
int imax(int i, int j)
Definition: TileLaserTimingTool.cxx:33
min
#define min(a, b)
Definition: cfImp.cxx:40
grepfile.ic
int ic
Definition: grepfile.py:33
create_dcsc_inputs_sqlite.arg
list arg
Definition: create_dcsc_inputs_sqlite.py:48
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
dumpNswErrorDb.maxStrip
tuple maxStrip
Definition: dumpNswErrorDb.py:27
CscFit.h
a
TList * a
Definition: liststreamerinfos.cxx:10
y
#define y
TH1F
Definition: rootspy.cxx:320
DeMoScan.first
bool first
Definition: DeMoScan.py:534
DEBUG
#define DEBUG
Definition: page_access.h:11
CscFit::cscfit
void cscfit(double *qstr, int &maxStrip, double &thr, double &da, int &ncl, double *sig, double *zpos, double &error) const
clustering by a Gaussing fit to the charge distribution 4 different fit are: First find the clusters ...
Definition: CscFit.cxx:30
CscFit::CscFit
CscFit(double sigma)
Definition: CscFit.cxx:20
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
merge.status
status
Definition: merge.py:17
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
WriteCellNoiseToCool.noise
noise
Definition: WriteCellNoiseToCool.py:380
Generated on Thu Jul 4 2024 21:13:14 for ATLAS Offline Software by   doxygen 1.8.18