All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
icarus::crt::CrtCal Class Reference

#include <CrtCal.h>

Public Member Functions

 CrtCal (const vector< TH1F * > *histos)
 
 ~CrtCal ()
 
void Cal ()
 
void PedCal ()
 
void GainCal ()
 
void GainFit (TH1F *h, size_t chan, float **statsarr, bool save)
 
void PedFit (TH1F *h, size_t chan, float *statsarr, bool save)
 
bool IsActive (TH1F *)
 
void ParsePedStats (const float *statarr, float &pedXsqr, float &ped, float &pedErr, float &pedNorm, float &pedNormErr, float &pedSigma, float &pedSigmaErr, short &pedNdf)
 
void ParseGainStats (float **statarr, float &gainXsqr, short &gainNdf, float &gain, float &gainErr, float &gainPed, float &gainPedErr, short &nPeak, float *peakXsqr, float *peakMean, float *peakMeanErr, float *peakNorm, float *peakNormErr, float *peakSigma, float *peakSigmaErr, short *peakNdf)
 
int FindThreshADC (TH1F *h)
 
float FindThreshPE (TH1F *h)
 
int FindNabove (TH1F *h, int thresh)
 
bool * GetActive () const
 
float * GetPed () const
 
float * GetPedErr () const
 
float * GetPedXsqr () const
 
short * GetPedNdf () const
 
float * GetPedNorm () const
 
float * GetPedNormErr () const
 
float * GetPedSigma () const
 
float * GetPedSigmaErr () const
 
float * GetGain () const
 
float * GetGainErr () const
 
float * GetGainXsqr () const
 
short * GetGainNdf () const
 
float * GetGainPed () const
 
float * GetGainPedErr () const
 
short * GetNpeak () const
 
int * GetThreshADC () const
 
float * GetThreshPE () const
 
int * GetNabove () const
 
float ** GetPeakNorm () const
 
float ** GetPeakNormErr () const
 
float ** GetPeakSigma () const
 
float ** GetPeakSigmaErr () const
 
float ** GetPeakMean () const
 
float ** GetPeakMeanErr () const
 
float ** GetPeakXsqr () const
 
short ** GetPeakNdf () const
 
long * GetChanStats () const
 
double * GetLangausWidth () const
 
double * GetLangausWidthErr () const
 
double * GetLangausLandauMP () const
 
double * GetLangausLandauMPErr () const
 
double * GetLangausArea () const
 
double * GetLangausAreaErr () const
 
double * GetLangausGaussSigma () const
 
double * GetLangausGaussSigmaErr () const
 
double * GetLangausXsqr () const
 
double * GetLangausNdf () const
 

Private Member Functions

void IndexToMacChan ()
 
TF1 * langaufit (TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
 
void langaus_fit (TH1F *h, double &lang_lan_wid, double &lang_lan_wid_err, double &lang_lan_mp, double &lang_lan_mp_err, double &lang_area, double &lang_area_err, double &lang_gauss_sigma, double &lang_gauss_sigma_err, double &lang_chisq, double &lang_ndf)
 

Private Attributes

const vector< TH1F * > * fHistos
 
map< uint8_t, uint8_t > fChanMap
 
bool fHasPedCal
 
bool fHasGainCal
 
bool fHasActive
 
bool fHasThresh
 
uint8_t fMac5
 
bool * fActive
 
float * fGain
 
float * fGainErr
 
float * fGainXsqr
 
short * fGainNdf
 
float * fGainPed
 
float * fGainPedErr
 
short * fNpeak
 
float ** fPeakNorm
 
float ** fPeakNormErr
 
float ** fPeakSigma
 
float ** fPeakSigmaErr
 
float ** fPeakMean
 
float ** fPeakMeanErr
 
float ** fPeakXsqr
 
short ** fPeakNdf
 
float * fPed
 
float * fPedErr
 
float * fPedXsqr
 
short * fPedNdf
 
float * fPedSigma
 
float * fPedSigmaErr
 
float * fPedNorm
 
float * fPedNormErr
 
int * fThreshADC
 
float * fThreshPE
 
int * fNabove
 
long * fChanStats
 
double * fLangausWidth
 
double * fLangausWidthErr
 
double * fLangausLandauMP
 
double * fLangausLandauMPErr
 
double * fLangausArea
 
double * fLangausAreaErr
 
double * fLangausGaussSigma
 
double * fLangausGaussSigmaErr
 
double * fLangausXsqr
 
double * fLangausNdf
 

Detailed Description

Definition at line 33 of file CrtCal.h.

Constructor & Destructor Documentation

CrtCal::CrtCal ( const vector< TH1F * > *  histos)

Definition at line 8 of file CrtCal.cc.

8  : fHistos(histos) {
9 
10  fMac5 = 0;
11  this->IndexToMacChan();
12 
13  fHasActive = false;
14  fHasThresh = false;
15  fHasPedCal = false;
16  fHasGainCal = false;
17 
18  fActive = new bool[32];
19 
20  fGain = new float[32];
21  fGainErr = new float[32];
22  fGainXsqr = new float[32];
23  fGainNdf = new short[32];
24  fGainPed = new float[32];
25  fGainPedErr = new float[32];
26  fNpeak = new short[32];
27 
28  fPed = new float[32];
29  fPedErr = new float[32];
30  fPedXsqr = new float[32];
31  fPedNdf = new short[32];
32  fPedSigma = new float[32];
33  fPedSigmaErr = new float[32];
34  fPedNorm = new float[32];
35  fPedNormErr = new float[32];
36  fThreshADC = new int[32];
37  fThreshPE = new float[32];
38  fNabove = new int[32];
39 
40  fChanStats = new long[32];
41 
42  fPeakNorm = new float*[32];
43  fPeakNormErr = new float*[32];
44  fPeakSigma = new float*[32];
45  fPeakSigmaErr = new float*[32];
46  fPeakMean = new float*[32];
47  fPeakMeanErr = new float*[32];
48  fPeakXsqr = new float*[32];
49  fPeakNdf = new short*[32];
50 
51  fLangausWidth = new double[32];
52  fLangausWidthErr = new double[32];
53  fLangausLandauMP = new double[32];
54  fLangausLandauMPErr = new double[32];
55  fLangausArea = new double[32];
56  fLangausAreaErr = new double[32];
57  fLangausGaussSigma = new double[32];
58  fLangausGaussSigmaErr = new double[32];
59  fLangausXsqr = new double[32];
60  fLangausNdf = new double[32];
61 
62  for(size_t i=0; i<32; i++){
63  fPeakNorm[i] = new float[5];
64  fPeakNormErr[i] = new float[5];
65  fPeakSigma[i] = new float[5];
66  fPeakSigmaErr[i] = new float[5];
67  fPeakMean[i] = new float[5];
68  fPeakMeanErr[i] = new float[5];
69  fPeakXsqr[i] = new float[5];
70  fPeakNdf[i] = new short[5];
71  }
72 
73 
74  for(size_t ch=0; ch<32; ch++){
75  fActive[ch] = false;
76 
77  fGain[ch] = FLT_MAX;
78  fGainErr[ch] = FLT_MAX;
79  fGainXsqr[ch] = FLT_MAX;
80  fGainNdf[ch] = SHRT_MAX;
81  fGainPed[ch] = FLT_MAX;
82  fGainPedErr[ch] = FLT_MAX;
83  fNpeak[ch] = SHRT_MAX;
84 
85  fPed[ch] = FLT_MAX;
86  fPedErr[ch] = FLT_MAX;
87  fPedXsqr[ch] = FLT_MAX;
88  fPedNdf[ch] = SHRT_MAX;
89  fPedSigma[ch] = FLT_MAX;
90  fPedSigmaErr[ch] = FLT_MAX;
91  fPedNorm[ch] = FLT_MAX;
92  fPedNormErr[ch] = FLT_MAX;
93  fThreshADC[ch] = INT_MAX;
94  fThreshPE[ch] = FLT_MAX;
95  fNabove[ch] = INT_MAX;
96  fChanStats[ch] = LONG_MAX;
97 
98  fLangausWidth[ch] = DBL_MAX;
99  fLangausWidthErr[ch] = DBL_MAX;
100  fLangausLandauMP[ch] = DBL_MAX;
101  fLangausLandauMPErr[ch] = DBL_MAX;
102  fLangausArea[ch] = DBL_MAX;
103  fLangausAreaErr[ch] = DBL_MAX;
104  fLangausGaussSigma[ch] = DBL_MAX;
105  fLangausGaussSigmaErr[ch] = DBL_MAX;
106  fLangausXsqr[ch] = DBL_MAX;
107  fLangausNdf[ch] = DBL_MAX;
108 
109 
110  for(size_t p=0; p<5; p++){
111  fPeakNorm[ch][p] = FLT_MAX;
112  fPeakNormErr[ch][p] = FLT_MAX;
113  fPeakSigma[ch][p] = FLT_MAX;
114  fPeakSigmaErr[ch][p] = FLT_MAX;
115  fPeakMean[ch][p] = FLT_MAX;
116  fPeakMeanErr[ch][p] = FLT_MAX;
117  fPeakXsqr[ch][p] = FLT_MAX;
118  fPeakNdf[ch][p] = SHRT_MAX;
119  }
120 
121  }
122 }
short * fPedNdf
Definition: CrtCal.h:134
float * fPedSigma
Definition: CrtCal.h:135
double * fLangausGaussSigma
Definition: CrtCal.h:150
float ** fPeakNorm
Definition: CrtCal.h:123
pdgs p
Definition: selectors.fcl:22
float * fPedNormErr
Definition: CrtCal.h:138
double * fLangausWidthErr
Definition: CrtCal.h:145
double * fLangausWidth
Definition: CrtCal.h:144
float ** fPeakSigmaErr
Definition: CrtCal.h:126
double * fLangausArea
Definition: CrtCal.h:148
double * fLangausLandauMP
Definition: CrtCal.h:146
float * fThreshPE
Definition: CrtCal.h:140
double * fLangausLandauMPErr
Definition: CrtCal.h:147
float * fGainErr
Definition: CrtCal.h:117
long * fChanStats
Definition: CrtCal.h:142
float * fPedSigmaErr
Definition: CrtCal.h:136
float ** fPeakMeanErr
Definition: CrtCal.h:128
float * fGainXsqr
Definition: CrtCal.h:118
float * fPedNorm
Definition: CrtCal.h:137
float ** fPeakSigma
Definition: CrtCal.h:125
short ** fPeakNdf
Definition: CrtCal.h:130
short * fGainNdf
Definition: CrtCal.h:119
float * fPedXsqr
Definition: CrtCal.h:133
double * fLangausAreaErr
Definition: CrtCal.h:149
const vector< TH1F * > * fHistos
Definition: CrtCal.h:104
float ** fPeakMean
Definition: CrtCal.h:127
double * fLangausNdf
Definition: CrtCal.h:153
float * fPedErr
Definition: CrtCal.h:132
float * fGainPedErr
Definition: CrtCal.h:121
void IndexToMacChan()
Definition: CrtCal.cc:181
double * fLangausGaussSigmaErr
Definition: CrtCal.h:151
float ** fPeakXsqr
Definition: CrtCal.h:129
float * fGainPed
Definition: CrtCal.h:120
float ** fPeakNormErr
Definition: CrtCal.h:124
short * fNpeak
Definition: CrtCal.h:122
double * fLangausXsqr
Definition: CrtCal.h:152
CrtCal::~CrtCal ( )

Definition at line 124 of file CrtCal.cc.

124  {
125 
126  delete[] fActive;
127 
128  delete[] fGain;
129  delete[] fGainErr;
130  delete[] fGainXsqr;
131  delete[] fGainNdf;
132  delete[] fGainPed;
133  delete[] fGainPedErr;
134  delete[] fNpeak;
135 
136  delete[] fPed;
137  delete[] fPedErr;
138  delete[] fPedXsqr;
139  delete[] fPedNdf;
140  delete[] fPedSigma;
141  delete[] fPedSigmaErr;
142  delete[] fPedNorm;
143  delete[] fPedNormErr;
144  delete[] fThreshADC;
145  delete[] fThreshPE;
146  delete[] fNabove;
147  delete[] fChanStats;
148 
149  for(size_t i=0; i<32; i++){
150  delete[] fPeakNorm[i];
151  delete[] fPeakNormErr[i];
152  delete[] fPeakSigma[i];
153  delete[] fPeakSigmaErr[i];
154  delete[] fPeakMean[i];
155  delete[] fPeakMeanErr[i];
156  delete[] fPeakXsqr[i];
157  delete[] fPeakNdf[i];
158  }
159 
160  delete[] fPeakNorm;
161  delete[] fPeakNormErr;
162  delete[] fPeakSigma;
163  delete[] fPeakSigmaErr;
164  delete[] fPeakMean;
165  delete[] fPeakMeanErr;
166  delete[] fPeakXsqr;
167  delete[] fPeakNdf;
168 
169  delete[] fLangausWidth;
170  delete[] fLangausWidthErr;
171  delete[] fLangausLandauMP;
172  delete[] fLangausArea;
173  delete[] fLangausGaussSigma;
174  delete[] fLangausLandauMPErr;
175  delete[] fLangausAreaErr;
176  delete[] fLangausGaussSigmaErr;
177  delete[] fLangausXsqr;
178  delete[] fLangausNdf;
179 }
short * fPedNdf
Definition: CrtCal.h:134
float * fPedSigma
Definition: CrtCal.h:135
double * fLangausGaussSigma
Definition: CrtCal.h:150
float ** fPeakNorm
Definition: CrtCal.h:123
float * fPedNormErr
Definition: CrtCal.h:138
double * fLangausWidthErr
Definition: CrtCal.h:145
double * fLangausWidth
Definition: CrtCal.h:144
float ** fPeakSigmaErr
Definition: CrtCal.h:126
double * fLangausArea
Definition: CrtCal.h:148
double * fLangausLandauMP
Definition: CrtCal.h:146
float * fThreshPE
Definition: CrtCal.h:140
double * fLangausLandauMPErr
Definition: CrtCal.h:147
float * fGainErr
Definition: CrtCal.h:117
long * fChanStats
Definition: CrtCal.h:142
float * fPedSigmaErr
Definition: CrtCal.h:136
float ** fPeakMeanErr
Definition: CrtCal.h:128
float * fGainXsqr
Definition: CrtCal.h:118
float * fPedNorm
Definition: CrtCal.h:137
float ** fPeakSigma
Definition: CrtCal.h:125
short ** fPeakNdf
Definition: CrtCal.h:130
short * fGainNdf
Definition: CrtCal.h:119
float * fPedXsqr
Definition: CrtCal.h:133
double * fLangausAreaErr
Definition: CrtCal.h:149
float ** fPeakMean
Definition: CrtCal.h:127
double * fLangausNdf
Definition: CrtCal.h:153
float * fPedErr
Definition: CrtCal.h:132
float * fGainPedErr
Definition: CrtCal.h:121
double * fLangausGaussSigmaErr
Definition: CrtCal.h:151
float ** fPeakXsqr
Definition: CrtCal.h:129
float * fGainPed
Definition: CrtCal.h:120
float ** fPeakNormErr
Definition: CrtCal.h:124
short * fNpeak
Definition: CrtCal.h:122
double * fLangausXsqr
Definition: CrtCal.h:152

Member Function Documentation

void CrtCal::Cal ( )

Definition at line 207 of file CrtCal.cc.

207  {
208 
209  size_t nactive = 0;
210  std::cout << "ActiveChannel scan..." << std::endl;
211  for(size_t i=0; i<fHistos->size(); i++){
212  TH1F* h1 = (TH1F*)fHistos->at(i)->Clone();
213  TH1F* h2 = (TH1F*)fHistos->at(i)->Clone();
214  TH1F* h3 = (TH1F*)fHistos->at(i)->Clone();
215  const size_t chan = fChanMap[i];
216  fActive[chan] = IsActive(h1);
217  delete h1;
218  if(fActive[chan]) {
219  nactive++;
220  fThreshADC[chan] = FindThreshADC(h2);
221  delete h2;
222  std::cout << "ch. " << chan << " ADC thresh: " << fThreshADC[chan] << std::endl;
223  fNabove[chan] = FindNabove(h3,fThreshADC[chan]);
224  delete h3;
225  }
226  }
227  std::cout << "found " << nactive << " active channels" << std::endl;
228 
229  std::cout << "PedCal..." << std::endl;
230  PedCal();
231  std::cout << "GainCal..." << std::endl;
232  GainCal();
233 
234  for(size_t ch=0; ch<32; ch++){
235  if(!fActive[ch]) continue;
236  fThreshPE[ch] = 1.0*(fThreshADC[ch]-fPed[ch])/fGain[ch];
237  }
238 }
map< uint8_t, uint8_t > fChanMap
Definition: CrtCal.h:107
float * fThreshPE
Definition: CrtCal.h:140
const vector< TH1F * > * fHistos
Definition: CrtCal.h:104
bool IsActive(TH1F *)
Definition: CrtCal.cc:240
int FindNabove(TH1F *h, int thresh)
Definition: CrtCal.cc:265
int FindThreshADC(TH1F *h)
Definition: CrtCal.cc:253
BEGIN_PROLOG could also be cout
int CrtCal::FindNabove ( TH1F *  h,
int  thresh 
)

Definition at line 265 of file CrtCal.cc.

265  {
266  return h->Integral(h->FindBin(thresh),h->GetNbinsX()+1);
267 }
while getopts h
int CrtCal::FindThreshADC ( TH1F *  h)

Definition at line 253 of file CrtCal.cc.

253  {
254  const size_t low = 300;
255  const size_t high = 1000;
256  int thresh=0;
257 
258  h->GetXaxis()->SetRangeUser(low,high);
259  thresh = (int)h->GetBinLowEdge(h->GetMinimumBin());
260  fHasThresh = true;
261 
262  return thresh;
263 }
standard_dbscan3dalg useful for diagnostics hits not in a line will not be clustered on on only for track like only for track like on on the smaller the less shower like tracks low
while getopts h
float icarus::crt::CrtCal::FindThreshPE ( TH1F *  h)
void CrtCal::GainCal ( )

Definition at line 293 of file CrtCal.cc.

293  {
294 
295  if(!fHasPedCal) {
296  std::cout << "need ped cal before gain cal!" << std::endl;
297  return;
298  }
299  if(fHasGainCal) return;
300 
301  float** statarr = new float*[15];
302  for(size_t i=0; i<15; i++){
303  statarr[i] = new float[5];
304  for(size_t j=0; j<5; j++){
305  statarr[i][j] = FLT_MAX;
306  }
307  }
308 
309  for(size_t i=0; i<fHistos->size(); i++){
310  const size_t chan = fChanMap[i];
311  if(!fActive[chan]) continue;
312  //std::cout << "channel " << chan << std::endl;
313  TH1F* h = (TH1F*)fHistos->at(i)->Clone();
314  fChanStats[chan]=h->GetEntries();
315  //std::cout << "call to GainFit()" << std::endl;
316  GainFit(h,chan,statarr,true);
317  //std::cout << "ParseGainStats()..." << std::endl;
318  ParseGainStats(statarr, fGainXsqr[chan], fGainNdf[chan], fGain[chan], fGainErr[chan],
319  fGainPed[chan], fGainPedErr[chan], fNpeak[chan], fPeakXsqr[chan],
320  fPeakMean[chan], fPeakMeanErr[chan], fPeakNorm[chan], fPeakNormErr[chan],
321  fPeakSigma[chan], fPeakSigmaErr[chan], fPeakNdf[chan]);
324  fLangausArea[chan],fLangausAreaErr[chan],
326  fLangausXsqr[chan],fLangausNdf[chan]);
327 
328  delete h;
329  }
330 
331  delete[] statarr;
332  fHasGainCal = true;
333  return;
334 }
map< uint8_t, uint8_t > fChanMap
Definition: CrtCal.h:107
double * fLangausGaussSigma
Definition: CrtCal.h:150
float ** fPeakNorm
Definition: CrtCal.h:123
void ParseGainStats(float **statarr, float &gainXsqr, short &gainNdf, float &gain, float &gainErr, float &gainPed, float &gainPedErr, short &nPeak, float *peakXsqr, float *peakMean, float *peakMeanErr, float *peakNorm, float *peakNormErr, float *peakSigma, float *peakSigmaErr, short *peakNdf)
Definition: CrtCal.cc:808
double * fLangausWidthErr
Definition: CrtCal.h:145
double * fLangausWidth
Definition: CrtCal.h:144
float ** fPeakSigmaErr
Definition: CrtCal.h:126
double * fLangausArea
Definition: CrtCal.h:148
double * fLangausLandauMP
Definition: CrtCal.h:146
void GainFit(TH1F *h, size_t chan, float **statsarr, bool save)
Definition: CrtCal.cc:453
double * fLangausLandauMPErr
Definition: CrtCal.h:147
while getopts h
float * fGainErr
Definition: CrtCal.h:117
long * fChanStats
Definition: CrtCal.h:142
float ** fPeakMeanErr
Definition: CrtCal.h:128
float * fGainXsqr
Definition: CrtCal.h:118
float ** fPeakSigma
Definition: CrtCal.h:125
short ** fPeakNdf
Definition: CrtCal.h:130
short * fGainNdf
Definition: CrtCal.h:119
double * fLangausAreaErr
Definition: CrtCal.h:149
const vector< TH1F * > * fHistos
Definition: CrtCal.h:104
float ** fPeakMean
Definition: CrtCal.h:127
double * fLangausNdf
Definition: CrtCal.h:153
float * fGainPedErr
Definition: CrtCal.h:121
double * fLangausGaussSigmaErr
Definition: CrtCal.h:151
float ** fPeakXsqr
Definition: CrtCal.h:129
float * fGainPed
Definition: CrtCal.h:120
void langaus_fit(TH1F *h, double &lang_lan_wid, double &lang_lan_wid_err, double &lang_lan_mp, double &lang_lan_mp_err, double &lang_area, double &lang_area_err, double &lang_gauss_sigma, double &lang_gauss_sigma_err, double &lang_chisq, double &lang_ndf)
Definition: CrtCal.cc:939
float ** fPeakNormErr
Definition: CrtCal.h:124
BEGIN_PROLOG could also be cout
short * fNpeak
Definition: CrtCal.h:122
double * fLangausXsqr
Definition: CrtCal.h:152
void CrtCal::GainFit ( TH1F *  h,
size_t  chan,
float **  statsarr,
bool  save 
)

Definition at line 453 of file CrtCal.cc.

453  {
454 
455  const float gainSeed = 55.0;
456  const Double_t tSpectrumSigma = 1.75;
457  const Double_t tSpectrumThreshold = 0.18;
458  const size_t nPeakMax = 5;
459 
460  h->Rebin(4);
461  h->GetXaxis()->SetRangeUser(350,700);
462  h->SetStats(kFALSE);
463 
464  TCanvas* cspec = new TCanvas();
465  h->Draw("e0hist");
466 
467  TSpectrum *s = new TSpectrum();
468  //args=source histo, sigma of searched peaks, options, threshold
469  int nPeak = s->Search(h,tSpectrumSigma,"",tSpectrumThreshold);
470  //std::cout << "found " << nPeak << " peaks to fit" << std::endl;
471 
472  //int ctr = 2;
473 
474  //std::cout << "declaring in function arrays" << std::endl;
475  float* peakXsqr = new float[5];
476  float* peakNdf = new float[5];
477  float* peakSigma = new float[5];
478  float* peakSigmaErr = new float[5];
479  float* peakNorm = new float[5];
480  float* peakNormErr = new float[5];
481  float* peakMean = new float[5];
482  float* peakMeanErr = new float[5];
483  /*float ** statarr = new float*[15];
484  for(size_t i=0; i<15; i++){
485  statarr[i] = new float[5];*/
486  for(size_t j=0; j<5; j++){
487  //statarr[i][j] = FLT_MAX;
488  peakXsqr[j] = FLT_MAX;
489  peakNdf[j] = FLT_MAX;
490  peakSigma[j] = FLT_MAX;
491  peakSigmaErr[j] = FLT_MAX;
492  peakNorm[j] = FLT_MAX;
493  peakNormErr[j] = FLT_MAX;
494  peakMean[j] = FLT_MAX;
495  peakMeanErr[j] = FLT_MAX;
496  }
497  //}
498 
499  //std::cout << "done." << std::endl;
500 
501  float *peds = GetPed();
502  float *pwidths = GetPedSigma();
503  //std::cout << "retreived ped info" << std::endl;
504 
505  Double_t *peaks = s->GetPositionX(); //candidate photopeaks ADC position
506 
507  //std::cout << "sorting peaks" << std::endl;
508  //Sort peaks
509  double temp = 0.;
510  int nchanges = 0;
511  do
512  {
513  nchanges=0;
514  for(int p=0;p<nPeak-1;p++)
515  if(peaks[p]>peaks[p+1])
516  {
517  temp=peaks[p];
518  peaks[p]=peaks[p+1];
519  peaks[p+1]=temp;
520  nchanges++;
521  }
522  } while( nchanges != 0 );
523 
524  //std::cout << "done." << std::endl;
525 
526  //ascending list of peak number(x) vs. ADC value(y) from TSpectrum
527  float x[10];
528  float y[10];//, ey[10];
529 
530  //same list with "bad" peaks excluded (what eventually goes into the gain fit)
531  float gx[11];
532  float gy[11], gey[11];
533  int peak_offset = round((peaks[0]-peds[chan])/gainSeed); //estimate peak number
534 
535  for (int j=0; j<nPeak&&j<10; j++)
536  {
537  if(peaks[j]<0) std::cout << "bad peak value: " << peaks[j] << std::endl;
538  x[j] = j+peak_offset;
539  y[j] = peaks[j];
540  }
541 
542  int gg = 1; //index of passing peak array
543  int nplow = 0; //no. peaks close to low hist edge
544  int nphigh = 0; //no. peaks close to high hist edge
545 
546  //first peak in passing array is pedestal
547  gx[0] = 0;
548  gy[0] = peds[chan];
549  gey[0] = pwidths[chan];
550 
551  //first and refit chi-squareds
552  double chisqr = 0.;//, chisqr0;
553  // int rnum = 0; //refit number
554 
555  //fit peaks about TSpectrum values
556  for (int g=0 ; g<nPeak&&g<(int)nPeakMax; g++)
557  {
558  //initial gaus fit to peak from TSpectrum
559  TF1 *gfit = new TF1("gfit","gaus",y[g]-20,y[g]+20);
560  gfit->SetParameter(0,h->GetBinContent(h->FindBin(y[g])));//200);
561  gfit->SetParameter(1,y[g]);
562  gfit->SetParameter(2,12);
563  gfit->SetParLimits(0,0,20000);
564  gfit->SetParLimits(1,y[g]-15,y[g]+15);
565  gfit->SetParLimits(2,8,40);
566  //std::cout << "fitting gaussian to peak detected at " << y[g] << " ADC..." << std::endl;
567 
568  h->Fit(gfit,"MQLR"); //fit peak
569  gfit->Draw("same");
570  //std::cout << "done" << std::endl;
571  //chisqr0 = gfit->GetChisquare()/gfit->GetNDF(); //get reduced chi-square from fit
572 
573  //check for peaks near the edges of the histogam
574  if( y[g]<h->GetBinLowEdge(1)+15) nplow++;
575  if( y[g]>h->GetBinLowEdge(h->GetNbinsX())-15) nphigh++;
576 
577  //ignore edge peaks and peaks with fit mean > 15 from peak
578  if( y[g]>h->GetBinLowEdge(1)+15 && y[g]<h->GetBinLowEdge(h->GetNbinsX())-15&&abs(y[g] - gfit->GetParameter(1)) < 15)
579  {
580  //skip false peaks (may need improvement - currently relies on init values)
581  if(g!=0 && g!=nPeak-1 //check it's not first or last peaks
582  && (y[g+1]-y[g]<gainSeed*0.7 || y[g]-y[g-1]<gainSeed*0.7))
583  { //check peak within 30% of expected gain w.r.t adj.
584  //std::cout << "determined peak " << g << " (" << y[g] << " ADC) is false peak. removing..." << std::endl;
585  for (int j=g; j<nPeak; j++) x[j] = x[j]-1;
586  }//overwrite current value, shifting higher values down by 1 index
587  //if not missed peak, could there have been a skipped peak?
588  else
589  {
590  if(g!=0&&g!=nPeak-1&&y[g+1]-y[g]<gainSeed*1.2&&y[g]-y[g-1]>gainSeed*1.5) //missed peak adjust
591  {
592  for (int j=g; j<nPeak; j++) x[j] = x[j]+1;
593  //std::cout << "missed peak detected. shifting peaks " << g << " and higher up by 1..." << std::endl;
594  }
595  //if last peak likely occuring after skipped peak
596  if(g==nPeak-1 && (y[g]-y[g-1])/gainSeed >1.5) {
597  x[g]+=(int)((y[g]-y[g-1])/gainSeed);
598  //std::cout << "missed peak(s) detected before last peak...shifting last peak" << std::endl;
599  }
600 
601  peakXsqr[gg] = gfit->GetChisquare();
602  peakNdf[gg] = gfit->GetNDF();
603  peakMean[gg] = gfit->GetParameter(1);
604  peakMeanErr[gg] = gfit->GetParError(1);
605  peakNorm[gg] = gfit->GetParameter(0);
606  peakNormErr[gg] = gfit->GetParError(0);
607  peakSigma[gg] = gfit->GetParameter(2);
608  peakSigmaErr[gg] = gfit->GetParError(2);
609 
610  gx[gg] = x[g];
611  gy[gg] = gfit->GetParameter(1);
612  gey[gg] = sqrt(gx[gg]+gey[0]*gey[0]);//gfit->GetParError(1);
613  gg++;
614  }//end else not missed but was there skip
615  }//end if not edge peak, not far from TSpectrum value
616  }//end for peaks
617 
618  //std::cout << "done looping over peaks" << std::endl;
619 
620  if (nplow>0) for (int i=1; i<gg; i++) gx[i] = gx[i]-(nplow-1);
621 
622  //graph of adc(y) vs. photo-peak number (x)
623  TGraphErrors* gr_mean = new TGraphErrors(gg,gx,gy,0,gey);
624 
625  //linear fit function
626  TF1 *fit = new TF1("fit","[0] + [1]*x",gx[0]-0.25,gx[gg-1]+0.25);
627 
628  //name, initialize gain fit parameters
629  fit->SetParName(1,"Gain");
630  fit->SetParName(0, "Pedestal");
631  fit->SetParameter(1,gainSeed);
632  fit->SetParameter(0,peds[chan]);
633  fit->SetParLimits(1,gainSeed-20,gainSeed+20);
634  fit->SetParLimits(0,peds[chan]*0.8,peds[chan]*1.2);
635 
636  //std::cout << "gain fit..." << std::endl;
637  //perform gain fit
638  gr_mean->Fit(fit, "QR");
639 
640  //std::cout << "check chi-square..." << std::endl;
641  //check if gain fit is bad according to chi-square
642  if (fit->GetChisquare()/fit->GetNDF()>5.0)
643  {
644  //std::cout << "gain fit X^2 too large...shifting all peaks by 1" << std::endl;
645  chisqr=fit->GetChisquare();
646  for(int i=1; i<gg; i++) gx[i]+=1;
647  gr_mean = new TGraphErrors(gg,gx,gy,0,gey);
648  fit->SetRange(gx[0]-0.25,gx[gg-1]+0.25);
649  gr_mean->Fit(fit, "QR");
650  if (fit->GetChisquare()<chisqr) chisqr=fit->GetChisquare();
651  else
652  {
653  for(int i=1; i<gg; i++) gx[i]-=2;
654  gr_mean = new TGraphErrors(gg,gx,gy,0,gey);
655  fit->SetRange(gx[0]-0.25,gx[gg-1]+0.25);
656  gr_mean->Fit(fit, "QR");
657  if (fit->GetChisquare()<chisqr) chisqr=fit->GetChisquare();
658  else
659  {
660  for(int i=1; i<gg; i++) gx[i]+=1;
661  gr_mean = new TGraphErrors(gg,gx,gy,0,gey);
662  fit->SetRange(gx[0]-0.25,gx[gg-1]+0.25);
663  gr_mean->Fit(fit, "QR");
664  }
665  }
666  }
667 
668  string grtitle = "Mac5 "+std::to_string(fMac5)+", Ch. "+std::to_string(chan)+" Gain Fit";
669 
670  //std::cout << "style and range opts" << std::endl;
671 
672  gr_mean->SetTitle(grtitle.c_str());
673  gr_mean->GetXaxis()->SetTitle("Peak #");
674  gr_mean->GetYaxis()->SetTitle("ADC value");
675  gr_mean->SetMarkerColor(4);
676  gr_mean->SetMarkerStyle(20);
677  gr_mean->SetFillColor(0);
678  //std::cout << "set range..." << std::endl;
679  gr_mean->GetXaxis()->SetRangeUser(gx[0]-0.5,gx[gg-1]+0.5);
680 
681  //std::cout << "done" << std::endl;
682  //if (fit->GetChisquare()/fit->GetNDF()>30.0) cout << "Poor X^2! mac5 " << mac << "ch. " << chan << endl;
683 
684  gStyle->SetOptStat(0100);
685  gStyle->SetOptFit(1111);
686 
687  //std::cout << "drawing..." << std::endl;
688 
689  TCanvas *c2 = new TCanvas();
690  c2->cd();
691  c2->SetGrid();
692 
693  gStyle->SetStatX(0.5);
694  gStyle->SetStatY(0.9);
695  gStyle->SetStatH(0.15);
696  gStyle->SetStatW(0.2);
697 
698  //std::cout << "draw gr_mean..." << std::endl;
699  gr_mean->Draw("ALP");
700  //std::cout << "draw fit" << std::endl;
701  fit->Draw("L,SAME");
702 
703  //std::cout << "filling statarr" << std::endl;
704 
705  statsarr[0][0] = fit->GetParameter(1); //gain
706  statsarr[1][0] = fit->GetParError(1); //gain error
707  statsarr[2][0] = fit->GetParameter(0); //pedestal mean
708  statsarr[3][0] = fit->GetParError(0); //pedestal mean error
709  statsarr[4][0] = fit->GetChisquare(); //X^2
710  statsarr[5][0] = fit->GetNDF(); //NDF
711  statsarr[6][0] = nPeak;
712  for(size_t i=0; i<5; i++) {
713  statsarr[7][i] = peakXsqr[i];
714  statsarr[8][i] = peakMean[i];
715  statsarr[9][i] = peakMeanErr[i];
716  statsarr[10][i] = peakNorm[i];
717  statsarr[11][i] = peakNormErr[i];
718  statsarr[12][i] = peakSigma[i];
719  statsarr[13][i] = peakSigmaErr[i];
720  statsarr[14][i] = peakNdf[i];
721  }
722 
723  delete[] peakXsqr;
724  delete[] peakMean;
725  delete[] peakMeanErr;
726  delete[] peakNorm;
727  delete[] peakNormErr;
728  delete[] peakSigma;
729  delete[] peakSigmaErr;
730  delete[] peakNdf;
731 
732  if ( save )
733  {
734  //std::cout << "saving plots..."<< std::endl;
735  int ctr=1, ctr2=1;
736  TString suff = ".png";
737  TString fname="./";
738  TString fname2 = fname;
739  fname+=fMac5;
740  fname2+=fMac5;
741  fname+="_ch";
742  fname2+="_ch";
743  if(chan<10) fname+="0";
744  if(chan<10) fname2+="0";
745  fname2+=chan;
746  fname2+="_peak_fit_spec";
747  fname+=to_string(chan)+"_fit-gain";
748  fname+="_1"+suff;
749  fname2+="_1"+suff;
750  while (!gSystem->AccessPathName(fname))
751  {
752  fname.Remove(fname.Length()-suff.Length(),fname.Length());
753  fname.Remove(fname.Length()-1,fname.Length());
754  ctr++;
755  fname+=to_string(ctr)+suff;
756  }
757 
758  while (!gSystem->AccessPathName(fname2))
759  {
760  fname2.Remove(fname2.Length()-suff.Length(),fname2.Length());
761  fname2.Remove(fname2.Length()-1,fname2.Length());
762  ctr2++;
763  fname2+=to_string(ctr2)+suff;
764  }
765 
766  //write image to file
767  TImage *img = TImage::Create();
768  img->FromPad(c2);
769  img->WriteImage(fname);
770 
771  //TImage *img2 = TImage::Create();
772  img->FromPad(cspec);
773  img->WriteImage(fname2);
774 
775  //deallocate memory
776  delete img;
777  delete c2;
778  delete cspec;
779  //delete img2;
780  delete gr_mean;
781  delete fit;
782  delete s;
783  }
784 
785  //std::cout << "done" << std::endl;
786 
787  return;// statarr;
788 
789 }
string fname
Definition: demo.py:5
process_name opflash particleana ie x
pdgs p
Definition: selectors.fcl:22
BEGIN_PROLOG g
float * GetPed() const
Definition: CrtCal.cc:972
while getopts h
T abs(T value)
process_name opflash particleana ie ie y
float * GetPedSigma() const
Definition: CrtCal.cc:1010
std::string to_string(WindowPattern const &pattern)
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
BEGIN_PROLOG could also be cout
bool * CrtCal::GetActive ( ) const

Definition at line 965 of file CrtCal.cc.

965  {
966  if(!fHasActive)
967  return nullptr;
968 
969  return fActive;
970 }
long * CrtCal::GetChanStats ( ) const

Definition at line 1083 of file CrtCal.cc.

1083  {
1084  if(!fHasGainCal)
1085  return nullptr;
1086 
1087  return fChanStats;
1088 }
long * fChanStats
Definition: CrtCal.h:142
float * CrtCal::GetGain ( ) const

Definition at line 1022 of file CrtCal.cc.

1022  {
1023  if(!fHasGainCal)
1024  return nullptr;
1025 
1026  return fGain;
1027 }
float * CrtCal::GetGainErr ( ) const

Definition at line 1028 of file CrtCal.cc.

1028  {
1029  if(!fHasGainCal)
1030  return nullptr;
1031 
1032  return fGainErr;
1033 }
float * fGainErr
Definition: CrtCal.h:117
short * CrtCal::GetGainNdf ( ) const

Definition at line 1040 of file CrtCal.cc.

1040  {
1041  if(!fHasGainCal)
1042  return nullptr;
1043 
1044  return fGainNdf;
1045 }
short * fGainNdf
Definition: CrtCal.h:119
float * CrtCal::GetGainPed ( ) const

Definition at line 1047 of file CrtCal.cc.

1047  {
1048  if(!fHasGainCal)
1049  return nullptr;
1050 
1051  return fGainPed;
1052 }
float * fGainPed
Definition: CrtCal.h:120
float * CrtCal::GetGainPedErr ( ) const

Definition at line 1053 of file CrtCal.cc.

1053  {
1054  if(!fHasGainCal)
1055  return nullptr;
1056 
1057  return fGainPedErr;
1058 }
float * fGainPedErr
Definition: CrtCal.h:121
float * CrtCal::GetGainXsqr ( ) const

Definition at line 1034 of file CrtCal.cc.

1034  {
1035  if(!fHasGainCal)
1036  return nullptr;
1037 
1038  return fGainXsqr;
1039 }
float * fGainXsqr
Definition: CrtCal.h:118
double * CrtCal::GetLangausArea ( ) const

Definition at line 1149 of file CrtCal.cc.

1149  {
1150  if(!fHasGainCal)
1151  return nullptr;
1152 
1153  return fLangausArea;
1154 }
double * fLangausArea
Definition: CrtCal.h:148
double * CrtCal::GetLangausAreaErr ( ) const

Definition at line 1173 of file CrtCal.cc.

1173  {
1174  if(!fHasGainCal)
1175  return nullptr;
1176 
1177  return fLangausAreaErr;
1178 }
double * fLangausAreaErr
Definition: CrtCal.h:149
double * CrtCal::GetLangausGaussSigma ( ) const

Definition at line 1155 of file CrtCal.cc.

1155  {
1156  if(!fHasGainCal)
1157  return nullptr;
1158 
1159  return fLangausGaussSigma;
1160 }
double * fLangausGaussSigma
Definition: CrtCal.h:150
double * CrtCal::GetLangausGaussSigmaErr ( ) const

Definition at line 1179 of file CrtCal.cc.

1179  {
1180  if(!fHasGainCal)
1181  return nullptr;
1182 
1183  return fLangausGaussSigmaErr;
1184 }
double * fLangausGaussSigmaErr
Definition: CrtCal.h:151
double * CrtCal::GetLangausLandauMP ( ) const

Definition at line 1143 of file CrtCal.cc.

1143  {
1144  if(!fHasGainCal)
1145  return nullptr;
1146 
1147  return fLangausLandauMP;
1148 }
double * fLangausLandauMP
Definition: CrtCal.h:146
double * CrtCal::GetLangausLandauMPErr ( ) const

Definition at line 1167 of file CrtCal.cc.

1167  {
1168  if(!fHasGainCal)
1169  return nullptr;
1170 
1171  return fLangausLandauMPErr;
1172 }
double * fLangausLandauMPErr
Definition: CrtCal.h:147
double * CrtCal::GetLangausNdf ( ) const

Definition at line 1191 of file CrtCal.cc.

1191  {
1192  if(!fHasGainCal)
1193  return nullptr;
1194 
1195  return fLangausNdf;
1196 }
double * fLangausNdf
Definition: CrtCal.h:153
double * CrtCal::GetLangausWidth ( ) const

Definition at line 1137 of file CrtCal.cc.

1137  {
1138  if(!fHasGainCal)
1139  return nullptr;
1140 
1141  return fLangausWidth;
1142 }
double * fLangausWidth
Definition: CrtCal.h:144
double * CrtCal::GetLangausWidthErr ( ) const

Definition at line 1161 of file CrtCal.cc.

1161  {
1162  if(!fHasGainCal)
1163  return nullptr;
1164 
1165  return fLangausWidthErr;
1166 }
double * fLangausWidthErr
Definition: CrtCal.h:145
double * CrtCal::GetLangausXsqr ( ) const

Definition at line 1185 of file CrtCal.cc.

1185  {
1186  if(!fHasGainCal)
1187  return nullptr;
1188 
1189  return fLangausXsqr;
1190 }
double * fLangausXsqr
Definition: CrtCal.h:152
int * CrtCal::GetNabove ( ) const

Definition at line 1077 of file CrtCal.cc.

1077  {
1078  if(!fHasThresh)
1079  return nullptr;
1080 
1081  return fNabove;
1082 }
short * CrtCal::GetNpeak ( ) const

Definition at line 1059 of file CrtCal.cc.

1059  {
1060  if(!fHasGainCal)
1061  return nullptr;
1062 
1063  return fNpeak;
1064 }
short * fNpeak
Definition: CrtCal.h:122
float ** CrtCal::GetPeakMean ( ) const

Definition at line 1113 of file CrtCal.cc.

1113  {
1114  if(!fHasGainCal)
1115  return nullptr;
1116 
1117  return fPeakMean;
1118 }
float ** fPeakMean
Definition: CrtCal.h:127
float ** CrtCal::GetPeakMeanErr ( ) const

Definition at line 1119 of file CrtCal.cc.

1119  {
1120  if(!fHasGainCal)
1121  return nullptr;
1122 
1123  return fPeakMeanErr;
1124 }
float ** fPeakMeanErr
Definition: CrtCal.h:128
short ** CrtCal::GetPeakNdf ( ) const

Definition at line 1131 of file CrtCal.cc.

1131  {
1132  if(!fHasGainCal)
1133  return nullptr;
1134 
1135  return fPeakNdf;
1136 }
short ** fPeakNdf
Definition: CrtCal.h:130
float ** CrtCal::GetPeakNorm ( ) const

Definition at line 1089 of file CrtCal.cc.

1089  {
1090  if(!fHasGainCal)
1091  return nullptr;
1092 
1093  return fPeakNorm;
1094 }
float ** fPeakNorm
Definition: CrtCal.h:123
float ** CrtCal::GetPeakNormErr ( ) const

Definition at line 1095 of file CrtCal.cc.

1095  {
1096  if(!fHasGainCal)
1097  return nullptr;
1098 
1099  return fPeakNormErr;
1100 }
float ** fPeakNormErr
Definition: CrtCal.h:124
float ** CrtCal::GetPeakSigma ( ) const

Definition at line 1101 of file CrtCal.cc.

1101  {
1102  if(!fHasGainCal)
1103  return nullptr;
1104 
1105  return fPeakSigma;
1106 }
float ** fPeakSigma
Definition: CrtCal.h:125
float ** CrtCal::GetPeakSigmaErr ( ) const

Definition at line 1107 of file CrtCal.cc.

1107  {
1108  if(!fHasGainCal)
1109  return nullptr;
1110 
1111  return fPeakSigmaErr;
1112 }
float ** fPeakSigmaErr
Definition: CrtCal.h:126
float ** CrtCal::GetPeakXsqr ( ) const

Definition at line 1125 of file CrtCal.cc.

1125  {
1126  if(!fHasGainCal)
1127  return nullptr;
1128 
1129  return fPeakXsqr;
1130 }
float ** fPeakXsqr
Definition: CrtCal.h:129
float * CrtCal::GetPed ( ) const

Definition at line 972 of file CrtCal.cc.

972  {
973  if(!fHasPedCal)
974  return nullptr;
975 
976  return fPed;
977 }
float * CrtCal::GetPedErr ( ) const

Definition at line 979 of file CrtCal.cc.

979  {
980  if(!fHasPedCal)
981  return nullptr;
982 
983  return fPedErr;
984 
985 }
float * fPedErr
Definition: CrtCal.h:132
short * CrtCal::GetPedNdf ( ) const

Definition at line 992 of file CrtCal.cc.

992  {
993  if(!fHasPedCal)
994  return nullptr;
995 
996  return fPedNdf;
997 }
short * fPedNdf
Definition: CrtCal.h:134
float * CrtCal::GetPedNorm ( ) const

Definition at line 998 of file CrtCal.cc.

998  {
999  if(!fHasPedCal)
1000  return nullptr;
1001 
1002  return fPedNorm;
1003 }
float * fPedNorm
Definition: CrtCal.h:137
float * CrtCal::GetPedNormErr ( ) const

Definition at line 1004 of file CrtCal.cc.

1004  {
1005  if(!fHasPedCal)
1006  return nullptr;
1007 
1008  return fPedNormErr;
1009 }
float * fPedNormErr
Definition: CrtCal.h:138
float * CrtCal::GetPedSigma ( ) const

Definition at line 1010 of file CrtCal.cc.

1010  {
1011  if(!fHasPedCal)
1012  return nullptr;
1013 
1014  return fPedSigma;
1015 }
float * fPedSigma
Definition: CrtCal.h:135
float * CrtCal::GetPedSigmaErr ( ) const

Definition at line 1016 of file CrtCal.cc.

1016  {
1017  if(!fHasPedCal)
1018  return nullptr;
1019 
1020  return fPedSigmaErr;
1021 }
float * fPedSigmaErr
Definition: CrtCal.h:136
float * CrtCal::GetPedXsqr ( ) const

Definition at line 986 of file CrtCal.cc.

986  {
987  if(!fHasPedCal)
988  return nullptr;
989 
990  return fPedXsqr;
991 }
float * fPedXsqr
Definition: CrtCal.h:133
int * CrtCal::GetThreshADC ( ) const

Definition at line 1065 of file CrtCal.cc.

1065  {
1066  if(!fHasThresh)
1067  return nullptr;
1068 
1069  return fThreshADC;
1070 }
float * CrtCal::GetThreshPE ( ) const

Definition at line 1071 of file CrtCal.cc.

1071  {
1072  if(!fHasGainCal)
1073  return nullptr;
1074 
1075  return fThreshPE;
1076 }
float * fThreshPE
Definition: CrtCal.h:140
void CrtCal::IndexToMacChan ( )
private

Definition at line 181 of file CrtCal.cc.

181  {
182 
183  if(fChanMap.size()!=0){
184  return;
185  }
186 
187  for(size_t i=0; i<fHistos->size(); i++){
188  //hname = hadc_mac_chan
189  const char* hname = fHistos->at(i)->GetName();
190  string name = "";
191  name+=hname;
192  if(i==0) {
193  const uint8_t nDigitMac = name.find("_",name.find("_")+1)-name.find("_")-1;
194  fMac5 = atoi(name.substr(name.find("_")+1,nDigitMac).c_str());
195  }
196  const uint8_t nDigitChan = name.size()-name.find("_",name.find("_")+1)-1;
197  uint8_t chan = atoi(name.substr(name.size()-nDigitChan,nDigitChan).c_str());
198  if(chan>31) std::cout << "ERROR in CrtCal::IndexToMacCHan: channel out of range (" << (short)chan << ")" << std::endl;
199  fChanMap[i] = chan;
200  }
201 
202  std::cout << "found " << fChanMap.size() << " histograms/channels" << std::endl;
203 
204  return;
205 }
map< uint8_t, uint8_t > fChanMap
Definition: CrtCal.h:107
const vector< TH1F * > * fHistos
Definition: CrtCal.h:104
then echo fcl name
BEGIN_PROLOG could also be cout
bool CrtCal::IsActive ( TH1F *  h)

Definition at line 240 of file CrtCal.cc.

240  {
241 
242  const size_t cutoff = 600;
243 
244  if(h->Integral(h->FindBin(cutoff),h->GetNbinsX())==0) {
245  return false;
246  }
247 
248  fHasActive = true;
249 
250  return true;
251 }
while getopts h
TF1 * CrtCal::langaufit ( TH1F *  his,
Double_t *  fitrange,
Double_t *  startvalues,
Double_t *  parlimitslo,
Double_t *  parlimitshi,
Double_t *  fitparams,
Double_t *  fiterrors,
Double_t *  ChiSqr,
Int_t *  NDF 
)
private

Definition at line 891 of file CrtCal.cc.

892  {
893  // Once again, here are the Landau * Gaussian parameters:
894  // par[0]=Width (scale) parameter of Landau density
895  // par[1]=Most Probable (MP, location) parameter of Landau density
896  // par[2]=Total area (integral -inf to inf, normalization constant)
897  // par[3]=Width (sigma) of convoluted Gaussian function
898  //
899  // Variables for langaufit call:
900  // his histogram to fit
901  // fitrange[2] lo and hi boundaries of fit range
902  // startvalues[4] reasonable start values for the fit
903  // parlimitslo[4] lower parameter limits
904  // parlimitshi[4] upper parameter limits
905  // fitparams[4] returns the final fit parameters
906  // fiterrors[4] returns the final fit errors
907  // ChiSqr returns the chi square
908  // NDF returns ndf
909 
910  Int_t i;
911  Char_t FunName[100];
912 
913  sprintf(FunName,"Fitfcn_%s",his->GetName());
914 
915 // TF1 *ffitold = (TF1*)ROOT->GetListOfFunctions()->FindObject(FunName);
916 // if (ffitold) delete ffitold;
917 
918  TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
919  ffit->SetParameters(startvalues);
920  ffit->SetParNames("Width","MP","Area","GSigma");
921 
922  for (i=0; i<4; i++) {
923  ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
924  }
925 
926  his->Fit(FunName,"RB0"); // fit within specified range, use ParLimits, do not plot
927 
928  ffit->GetParameters(fitparams); // obtain fit parameters
929  for (i=0; i<4; i++) {
930  fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors
931  }
932  ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2
933  NDF[0] = ffit->GetNDF(); // obtain ndf
934 
935  return (ffit); // return fit function
936 
937  }
Double_t langaufun(Double_t *x, Double_t *par)
Definition: CrtCal.cc:835
void CrtCal::langaus_fit ( TH1F *  h,
double &  lang_lan_wid,
double &  lang_lan_wid_err,
double &  lang_lan_mp,
double &  lang_lan_mp_err,
double &  lang_area,
double &  lang_area_err,
double &  lang_gauss_sigma,
double &  lang_gauss_sigma_err,
double &  lang_chisq,
double &  lang_ndf 
)
private

Definition at line 939 of file CrtCal.cc.

939  {
940 
941  //set up bounds, initial guesses for fit
942  Double_t my_fitrange[2] = { 800,4000 };
943  Double_t my_startvalues[4] = { 3, 200, 10000, 150 };
944  Double_t my_parlimitslo[4] = { 0, 100, 1000, 100 };
945  Double_t my_parlimitshi[4] = { 1000, 10000, 10000000, 1000 };
946 
947  //set up receptacles for relevant fit information
948  Double_t my_fitparams[4], my_fiterrors[4], my_ChiSqr;
949  Int_t my_NDF;
950 
951  TF1 *my_langaus = langaufit(h, my_fitrange, my_startvalues, my_parlimitslo, my_parlimitshi, my_fitparams, my_fiterrors, &my_ChiSqr, &my_NDF);
952 
953  lang_lan_wid = my_fitparams[0]; lang_lan_wid_err = my_fiterrors[0];
954  lang_lan_mp = my_fitparams[1]; lang_lan_mp_err = my_fiterrors[1];
955  lang_area = my_fitparams[2]; lang_area_err = my_fiterrors[2];
956  lang_gauss_sigma = my_fitparams[3]; lang_gauss_sigma_err = my_fiterrors[3];
957  lang_chisq = my_ChiSqr; lang_ndf = my_NDF;
958 
959 
960  delete my_langaus;
961 }
TF1 * langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
Definition: CrtCal.cc:891
while getopts h
void CrtCal::ParseGainStats ( float **  statarr,
float &  gainXsqr,
short &  gainNdf,
float &  gain,
float &  gainErr,
float &  gainPed,
float &  gainPedErr,
short &  nPeak,
float *  peakXsqr,
float *  peakMean,
float *  peakMeanErr,
float *  peakNorm,
float *  peakNormErr,
float *  peakSigma,
float *  peakSigmaErr,
short *  peakNdf 
)

Definition at line 808 of file CrtCal.cc.

811  {
812 
813  gain = statarr[0][0];
814  gainErr = statarr[1][0];
815  gainPed = statarr[2][0];
816  gainPedErr = statarr[3][0];
817  gainXsqr = statarr[4][0];
818  gainNdf = (short)statarr[5][0];
819  nPeak = statarr[6][0];
820  for(size_t i=0; i<5; i++){
821  //std::cout << "i="<< i << std::endl;
822  peakXsqr[i] = statarr[7][i];
823  peakMean[i] = statarr[8][i];
824  peakMeanErr[i] = statarr[9][i];
825  peakNorm[i] = statarr[10][i];
826  peakNormErr[i] = statarr[11][i];
827  peakSigma[i] = statarr[12][i];
828  peakSigmaErr[i] = statarr[13][i];
829  peakNdf[i] = (short)statarr[14][i];
830  }
831 
832  return;
833 }
void CrtCal::ParsePedStats ( const float *  statarr,
float &  pedXsqr,
float &  ped,
float &  pedErr,
float &  pedNorm,
float &  pedNormErr,
float &  pedSigma,
float &  pedSigmaErr,
short &  pedNdf 
)

Definition at line 792 of file CrtCal.cc.

794  {
795 
796  pedNorm = statarr[0]; //const
797  pedNormErr = statarr[1]; //const err
798  ped = statarr[2]; //mean
799  pedErr = statarr[3]; //mean err
800  pedSigma = statarr[4]; //sigma
801  pedSigmaErr = statarr[5]; //sigma err
802  pedXsqr = statarr[6];
803  pedNdf = (short)statarr[7];
804 
805  return;
806 }
void CrtCal::PedCal ( )

Definition at line 269 of file CrtCal.cc.

269  {
270 
271  if(fHasPedCal) return;
272 
273  float* statarr = new float[8];
274  for(size_t i=0; i<8; i++){
275  statarr[i] = FLT_MAX;
276  }
277 
278  for(size_t i=0; i<fHistos->size(); i++){
279  const size_t chan = fChanMap[i];
280  TH1F* h = (TH1F*)fHistos->at(i)->Clone();
281  PedFit(h,chan,statarr,true);
282  ParsePedStats(statarr, fPed[chan], fPedErr[chan], fPedNorm[chan], fPedNormErr[chan],
283  fPedSigma[chan], fPedSigmaErr[chan], fPedXsqr[chan], fPedNdf[chan] );
284  delete h;
285  }
286 
287  delete[] statarr;
288  fHasPedCal = true;
289 
290  return;
291 }
short * fPedNdf
Definition: CrtCal.h:134
float * fPedSigma
Definition: CrtCal.h:135
map< uint8_t, uint8_t > fChanMap
Definition: CrtCal.h:107
float * fPedNormErr
Definition: CrtCal.h:138
void PedFit(TH1F *h, size_t chan, float *statsarr, bool save)
Definition: CrtCal.cc:336
while getopts h
float * fPedSigmaErr
Definition: CrtCal.h:136
float * fPedNorm
Definition: CrtCal.h:137
float * fPedXsqr
Definition: CrtCal.h:133
const vector< TH1F * > * fHistos
Definition: CrtCal.h:104
void ParsePedStats(const float *statarr, float &pedXsqr, float &ped, float &pedErr, float &pedNorm, float &pedNormErr, float &pedSigma, float &pedSigmaErr, short &pedNdf)
Definition: CrtCal.cc:792
float * fPedErr
Definition: CrtCal.h:132
void CrtCal::PedFit ( TH1F *  h,
size_t  chan,
float *  statsarr,
bool  save = false 
)

Definition at line 336 of file CrtCal.cc.

336  {
337 
338  const int low = 1, high = 300;
339 
340  h->GetXaxis()->SetRangeUser(low,high);
341 
342  const int max_bin = h->GetMaximumBin();
343  const int max_val = h->GetMaximum();
344  const int max_adc = h->GetBinLowEdge(max_bin);
345 
346  //Define fitting function:
347  TF1 *gausfit = new TF1("gausfit","[0]*exp(-0.5*((x-[1])/[2])^2)",max_adc-12,max_adc+12);
348 
349  //Set parameter names:
350  gausfit->SetParName(0,"Constant");
351  gausfit->SetParName(1,"Peak value");
352  gausfit->SetParName(2,"sigma");
353 
354  //Initial guesses for the parameters:
355  gausfit->SetParameter(0,max_val);
356  gausfit->SetParLimits(0,0.5*max_val,1000*max_val);
357  gausfit->SetParameter(1,max_adc);
358  gausfit->SetParLimits(1,max_adc-20,max_adc+20);
359  gausfit->SetParameter(2,50.0);//h->GetStdDev());
360  gausfit->SetParLimits(2,1,50);
361 
362  TCanvas *c = new TCanvas();
363  c->SetGrid();
364 
365  gStyle->Reset("Modern");
366  if(h->GetMean()>(high-low)/2)
367  gStyle->SetStatX(0.4);
368  else
369  gStyle->SetStatX(0.9);
370 
371  gStyle->SetStatY(0.89);
372  gStyle->SetStatH(0.45);
373  gStyle->SetStatW(0.28);
374  h->Draw("e0samehist");
375  //h->Draw("E0same");//hist");
376  h->Fit(gausfit,"RMQ");
377 
378  float csq = gausfit->GetChisquare(); //chi-squared
379  float ndf = gausfit->GetNDF(); //NDF
380  float rcsq = csq/ndf; //reduced chi-squared
381 
382  if(rcsq>10.0)
383  {
384  //cout << "X^2/NDF > 10...recursive refit..." << endl;
385  gausfit->SetRange(gausfit->GetParameter(1)-10,gausfit->GetParameter(1)+10);
386  gausfit->SetParameter(0,gausfit->GetParameter(0));
387  gausfit->SetParLimits(0,gausfit->GetParameter(0)-50,gausfit->GetParameter(0)+50);
388  gausfit->SetParameter(1,gausfit->GetParameter(1));
389  gausfit->SetParLimits(1,gausfit->GetParameter(1)-10,gausfit->GetParameter(1)+10);
390  gausfit->SetParameter(2,gausfit->GetParameter(2));
391  gausfit->SetParLimits(2,gausfit->GetParameter(2)-10,gausfit->GetParameter(2)+10);
392  h->Fit(gausfit,"RMQ");
393  }
394 
395  gausfit->Draw("same");
396 
397  //float* statsarr = new float[8];
398  statsarr[0]=gausfit->GetParameter(0); //const
399  statsarr[1]=gausfit->GetParError(0); //const err
400  statsarr[2]=gausfit->GetParameter(1); //mean
401  statsarr[3]=gausfit->GetParError(1); //mean err
402  statsarr[4]=gausfit->GetParameter(2); //sigma
403  statsarr[5]=gausfit->GetParError(2); //sigma err
404  statsarr[6]=gausfit->GetChisquare();
405  statsarr[7]=gausfit->GetNDF();
406 
407  //std::cout << " ped value = " << statarr[2] << std::endl;
408 
409  //Adding chi^2/ndf and fit parameters to stat box:
410  gStyle->SetOptStat(1111);
411  gStyle->SetOptFit(111);
412  c->Update();
413 
414  csq = gausfit->GetChisquare();
415  ndf = gausfit->GetNDF();
416  rcsq = csq/ndf;
417  if (rcsq>200.0) std::cout << "warning: possibly bad ped fit mac5 " << fMac5 << ", ch. "
418  << chan << " X^2/NDF=" << rcsq << " ADC" << std::endl;
419  //save plot to file if desired
420  if ( save )
421  {
422  int ctr = 1;
423  TString suff=".png";
424  TString fname = "./";//calPlotDir+"/";
425  fname+=fMac5;
426  fname+="_";
427  fname+="chan"; //FormatMacString(mac)+ "_ch";
428  if(chan<10) fname+="0";
429  fname+=to_string(chan)+"_pedfit_"+to_string(ctr)+suff;
430 
431  while (!gSystem->AccessPathName(fname))
432  {
433  fname.Remove(fname.Length()-suff.Length(),fname.Length());
434  fname.Remove(fname.Length()-1,fname.Length());
435  ctr++;
436  fname+=to_string(ctr)+suff;
437  }
438 
439  //write image to file
440  TImage *img = TImage::Create();
441  img->FromPad(c);
442  img->WriteImage(fname);
443 
444  delete img;
445  delete gausfit;
446  delete c;
447  }//end save opt
448 
449  return;// statarr;
450 
451 }
string fname
Definition: demo.py:5
standard_dbscan3dalg useful for diagnostics hits not in a line will not be clustered on on only for track like only for track like on on the smaller the less shower like tracks low
while getopts h
std::string to_string(WindowPattern const &pattern)
BEGIN_PROLOG could also be cout

Member Data Documentation

bool* icarus::crt::CrtCal::fActive
private

Definition at line 115 of file CrtCal.h.

map<uint8_t,uint8_t> icarus::crt::CrtCal::fChanMap
private

Definition at line 107 of file CrtCal.h.

long* icarus::crt::CrtCal::fChanStats
private

Definition at line 142 of file CrtCal.h.

float* icarus::crt::CrtCal::fGain
private

Definition at line 116 of file CrtCal.h.

float* icarus::crt::CrtCal::fGainErr
private

Definition at line 117 of file CrtCal.h.

short* icarus::crt::CrtCal::fGainNdf
private

Definition at line 119 of file CrtCal.h.

float* icarus::crt::CrtCal::fGainPed
private

Definition at line 120 of file CrtCal.h.

float* icarus::crt::CrtCal::fGainPedErr
private

Definition at line 121 of file CrtCal.h.

float* icarus::crt::CrtCal::fGainXsqr
private

Definition at line 118 of file CrtCal.h.

bool icarus::crt::CrtCal::fHasActive
private

Definition at line 111 of file CrtCal.h.

bool icarus::crt::CrtCal::fHasGainCal
private

Definition at line 110 of file CrtCal.h.

bool icarus::crt::CrtCal::fHasPedCal
private

Definition at line 109 of file CrtCal.h.

bool icarus::crt::CrtCal::fHasThresh
private

Definition at line 112 of file CrtCal.h.

const vector<TH1F*>* icarus::crt::CrtCal::fHistos
private

Definition at line 104 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausArea
private

Definition at line 148 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausAreaErr
private

Definition at line 149 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausGaussSigma
private

Definition at line 150 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausGaussSigmaErr
private

Definition at line 151 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausLandauMP
private

Definition at line 146 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausLandauMPErr
private

Definition at line 147 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausNdf
private

Definition at line 153 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausWidth
private

Definition at line 144 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausWidthErr
private

Definition at line 145 of file CrtCal.h.

double* icarus::crt::CrtCal::fLangausXsqr
private

Definition at line 152 of file CrtCal.h.

uint8_t icarus::crt::CrtCal::fMac5
private

Definition at line 114 of file CrtCal.h.

int* icarus::crt::CrtCal::fNabove
private

Definition at line 141 of file CrtCal.h.

short* icarus::crt::CrtCal::fNpeak
private

Definition at line 122 of file CrtCal.h.

float** icarus::crt::CrtCal::fPeakMean
private

Definition at line 127 of file CrtCal.h.

float** icarus::crt::CrtCal::fPeakMeanErr
private

Definition at line 128 of file CrtCal.h.

short** icarus::crt::CrtCal::fPeakNdf
private

Definition at line 130 of file CrtCal.h.

float** icarus::crt::CrtCal::fPeakNorm
private

Definition at line 123 of file CrtCal.h.

float** icarus::crt::CrtCal::fPeakNormErr
private

Definition at line 124 of file CrtCal.h.

float** icarus::crt::CrtCal::fPeakSigma
private

Definition at line 125 of file CrtCal.h.

float** icarus::crt::CrtCal::fPeakSigmaErr
private

Definition at line 126 of file CrtCal.h.

float** icarus::crt::CrtCal::fPeakXsqr
private

Definition at line 129 of file CrtCal.h.

float* icarus::crt::CrtCal::fPed
private

Definition at line 131 of file CrtCal.h.

float* icarus::crt::CrtCal::fPedErr
private

Definition at line 132 of file CrtCal.h.

short* icarus::crt::CrtCal::fPedNdf
private

Definition at line 134 of file CrtCal.h.

float* icarus::crt::CrtCal::fPedNorm
private

Definition at line 137 of file CrtCal.h.

float* icarus::crt::CrtCal::fPedNormErr
private

Definition at line 138 of file CrtCal.h.

float* icarus::crt::CrtCal::fPedSigma
private

Definition at line 135 of file CrtCal.h.

float* icarus::crt::CrtCal::fPedSigmaErr
private

Definition at line 136 of file CrtCal.h.

float* icarus::crt::CrtCal::fPedXsqr
private

Definition at line 133 of file CrtCal.h.

int* icarus::crt::CrtCal::fThreshADC
private

Definition at line 139 of file CrtCal.h.

float* icarus::crt::CrtCal::fThreshPE
private

Definition at line 140 of file CrtCal.h.


The documentation for this class was generated from the following files: