ATLAS Offline Software
Public Member Functions | Static Private Member Functions | Private Attributes | Friends | List of all members
CscFit Class Reference

#include <CscFit.h>

Collaboration diagram for CscFit:

Public Member Functions

 CscFit (double sigma)
 
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 for each cluster fit to the left and to the right of the peak then fit the whole distribution with one Gaussian Finally fit the whold distribution with 2 Gaussians If the fit does not converge or is "bad" then take the weighted average More...
 

Static Private Member Functions

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 More...
 

Private Attributes

double m_sigma
 

Friends

Double_t f1gauss (const Double_t *x, const Double_t *par)
 one Gaussing fitting method More...
 
Double_t f2gauss (const Double_t *x, const Double_t *par)
 two-Gaussian fitting method More...
 

Detailed Description

Definition at line 20 of file CscFit.h.

Constructor & Destructor Documentation

◆ CscFit()

CscFit::CscFit ( double  sigma)

Definition at line 20 of file CscFit.cxx.

20 { m_sigma = sigma; }

Member Function Documentation

◆ cscfit()

void CscFit::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 for each cluster fit to the left and to the right of the peak then fit the whole distribution with one Gaussian Finally fit the whold distribution with 2 Gaussians If the fit does not converge or is "bad" then take the weighted average

A Gaussin fit to the left of the strip with the maximum charge in the cluster

A Gaussin fit to the right of the strip with the maximum charge in the cluster

One Gaussin fit to the cluster charge distribution

Two-Gaussin fit to the cluster charge distribution

the error on the weighted mean!!

make a selection for the best position finding the error on the position is the error on the weight mean

Definition at line 30 of file CscFit.cxx.

30  {
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 }

◆ icmax()

int CscFit::icmax ( const double *  qstr,
const int &  i1,
const int &  i2 
)
staticprivate

a method to find a the index of the strip with the highest charge

Definition at line 334 of file CscFit.cxx.

334  {
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 }

Friends And Related Function Documentation

◆ f1gauss

Double_t f1gauss ( const Double_t *  x,
const Double_t *  par 
)
friend

one Gaussing fitting method

Definition at line 349 of file CscFit.cxx.

349  {
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 }

◆ f2gauss

Double_t f2gauss ( const Double_t *  x,
const Double_t *  par 
)
friend

two-Gaussian fitting method

Definition at line 357 of file CscFit.cxx.

357  {
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 }

Member Data Documentation

◆ m_sigma

double CscFit::m_sigma
private

Definition at line 35 of file CscFit.h.

The documentation for this class was generated from the following files:
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
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
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
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
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
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 Mon Jul 8 2024 21:42:30 for ATLAS Offline Software by   doxygen 1.8.18