source: trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityRelTimeCam.cc@ 5974

Last change on this file since 5974 was 5812, checked in by gaug, 21 years ago
*** empty log message ***
File size: 9.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Markus Gaug 11/2003 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24/////////////////////////////////////////////////////////////////////////////
25//
26// MCalibrationIntensityRelTimeCam
27//
28// Storage container for intensity charge calibration results.
29//
30// Individual MCalibrationRelTimeCam's can be retrieved with:
31// - GetCam() yielding the current cam.
32// - GetCam("name") yielding the current camera with name "name".
33// - GetCam(i) yielding the i-th camera.
34//
35// See also: MCalibrationIntensityCam, MCalibrationRelTimeCam,
36// MCalibrationRelTimePix, MCalibrationRelTimeCalc, MCalibrationQECam
37// MHCalibrationRelTimePix, MHCalibrationRelTimeCam
38//
39/////////////////////////////////////////////////////////////////////////////
40#include "MCalibrationIntensityRelTimeCam.h"
41#include "MCalibrationRelTimeCam.h"
42#include "MCalibrationChargeCam.h"
43#include "MCalibrationRelTimePix.h"
44#include "MCalibrationChargePix.h"
45
46#include "MGeomCam.h"
47#include "MGeomPix.h"
48
49#include "MHCamera.h"
50
51#include "MLogManip.h"
52
53#include <TOrdCollection.h>
54#include <TGraphErrors.h>
55#include <TH1D.h>
56#include <TF1.h>
57
58ClassImp(MCalibrationIntensityRelTimeCam);
59
60using namespace std;
61
62// --------------------------------------------------------------------------
63//
64// Default constructor.
65//
66MCalibrationIntensityRelTimeCam::MCalibrationIntensityRelTimeCam(const char *name, const char *title)
67{
68
69 fName = name ? name : "MCalibrationIntensityRelTimeCam";
70 fTitle = title ? title : "Results of the Intensity Calibration";
71
72 InitSize(1);
73}
74
75// -------------------------------------------------------------------
76//
77// Add MCalibrationRelTimeCam's in the ranges from - to.
78//
79void MCalibrationIntensityRelTimeCam::Add(const UInt_t from, const UInt_t to)
80{
81 for (UInt_t i=from; i<to; i++)
82 fCams->AddAt(new MCalibrationRelTimeCam,i);
83}
84
85TGraphErrors *MCalibrationIntensityRelTimeCam::GetTimeResVsCharge( const UInt_t pixid, const MCalibrationIntensityChargeCam &chargecam,
86 const MCalibrationCam::PulserColor_t col)
87{
88
89 if (chargecam.GetSize() != GetSize())
90 {
91 *fLog << err << GetDescriptor() << ": Size mismatch between MCalibrationIntensityRelTimeCam "
92 << "and MCalibrationIntensityChargeCam. " << endl;
93 return NULL;
94 }
95
96 Int_t size = CountNumEntries(col);
97
98 if (size == 0)
99 return NULL;
100
101 if (size != chargecam.CountNumEntries(col))
102 {
103 *fLog << err << GetDescriptor() << ": Size mismatch in colour between MCalibrationIntensityRelTimeCam "
104 << "and MCalibrationIntensityChargeCam. " << endl;
105 return NULL;
106 }
107
108 const Int_t nvalid = chargecam.CountNumValidEntries(pixid,col);
109
110 if (nvalid == 0)
111 {
112 *fLog << err << GetDescriptor() << ": Only un-calibrated events in pixel: " << pixid << endl;
113 return NULL;
114 }
115
116 TArrayF res(nvalid);
117 TArrayF reserr(nvalid);
118 TArrayF sig(nvalid);
119 TArrayF sigerr(nvalid);
120
121 const Float_t sqrt2 = 1.414;
122 const Float_t fadc2ns = 3.333;
123
124 Int_t cnt = 0;
125
126 for (Int_t i=0;i<GetSize();i++)
127 {
128 //
129 // Get the calibration cam from the intensity cam
130 //
131 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)chargecam.GetCam(i);
132 MCalibrationRelTimeCam *relcam = (MCalibrationRelTimeCam*)GetCam(i);
133
134 if (col != MCalibrationCam::kNONE)
135 if (relcam->GetPulserColor() != col)
136 continue;
137 //
138 // Get the calibration pix from the calibration cam
139 //
140 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
141 MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[pixid];
142 //
143 // Don't use bad pixels
144 //
145 if (!pix.IsFFactorMethodValid())
146 continue;
147 //
148 res[cnt] = relpix.GetTimePrecision() / sqrt2 * fadc2ns;
149 reserr[cnt] = relpix.GetTimePrecisionErr() / sqrt2 * fadc2ns;
150 //
151 sig [cnt] = pix.GetPheFFactorMethod();
152 sigerr[cnt] = pix.GetPheFFactorMethodErr();
153 cnt++;
154 }
155
156 TGraphErrors *gr = new TGraphErrors(nvalid,
157 sig.GetArray(),res.GetArray(),
158 sigerr.GetArray(),reserr.GetArray());
159 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
160 gr->GetXaxis()->SetTitle("<Photo-electrons> [1]");
161 gr->GetYaxis()->SetTitle("Time Resolution [ns]");
162 return gr;
163}
164
165
166TGraphErrors *MCalibrationIntensityRelTimeCam::GetTimeResVsChargePerArea( const Int_t aidx,const MCalibrationIntensityChargeCam &chargecam, const MGeomCam &geom, const MCalibrationCam::PulserColor_t col)
167{
168
169 Int_t size = CountNumEntries(col);
170
171 if (size == 0)
172 return NULL;
173
174 if (size != chargecam.CountNumEntries(col))
175 {
176 *fLog << err << GetDescriptor() << ": Size mismatch in colour between MCalibrationIntensityRelTimeCam "
177 << "and MCalibrationIntensityChargeCam. " << endl;
178 return NULL;
179 }
180
181 const Float_t sqrt2 = 1.414;
182 const Float_t fadc2ns = 3.333;
183 const Float_t norm = fadc2ns / sqrt2;
184
185 TArrayF res(size);
186 TArrayF reserr(size);
187 TArrayF sig(size);
188 TArrayF sigerr(size);
189
190 Int_t cnt = 0;
191
192 TH1D *h = 0;
193
194 for (Int_t i=0;i<GetSize();i++)
195 {
196 //
197 // Get the calibration cam from the intensity cam
198 //
199 MCalibrationRelTimeCam *relcam = (MCalibrationRelTimeCam*)GetCam(i);
200 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)chargecam.GetCam(i);
201
202 if (relcam->GetPulserColor() != col)
203 continue;
204
205 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
206 const Float_t phe = apix.GetPheFFactorMethod();
207 const Float_t pheerr = apix.GetPheFFactorMethodErr();
208
209 sig[cnt] = phe;
210 sigerr[cnt] = pheerr;
211
212 Double_t resol = 0.;
213 Double_t resol2 = 0.;
214 Double_t var = 0.;
215 Int_t num = 0;
216
217 MHCamera camres(geom,"CamRes","Time Resolution;Time Resolution [ns];channels");
218 //
219 // Get the area calibration pix from the calibration cam
220 //
221 for (Int_t j=0; j<cam->GetSize(); j++)
222 {
223 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[j];
224 const MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[j];
225 //
226 // Don't use bad pixels
227 //
228 if (!pix.IsFFactorMethodValid())
229 continue;
230 //
231 //
232 if (aidx != geom[j].GetAidx())
233 continue;
234
235 const Float_t pres = relpix.GetTimePrecision();
236
237 resol += pres;
238 resol2 += pres;
239 num++;
240
241 camres.Fill(j,pres);
242 camres.SetUsed(j);
243 }
244
245 if (num > 1)
246 {
247 resol /= num;
248 var = (resol2 - resol*resol*num) / (num-1);
249
250 res[cnt] = resol;
251 if (var > 0.)
252 reserr[cnt] = TMath::Sqrt(var);
253 else
254 reserr[cnt] = 0.;
255
256 //
257 // Make also a Gauss-fit to the distributions. The RMS can be determined by
258 // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
259 //
260 h = camres.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
261 h->SetDirectory(NULL);
262 h->Fit("gaus","QL");
263 TF1 *fit = h->GetFunction("gaus");
264
265 Float_t ci2 = fit->GetChisquare();
266 Float_t sigma = fit->GetParameter(2);
267
268 if (ci2 > 500. || sigma > reserr[cnt])
269 {
270 h->Fit("gaus","QLM");
271 fit = h->GetFunction("gaus");
272
273 ci2 = fit->GetChisquare();
274 sigma = fit->GetParameter(2);
275 }
276
277 const Float_t mean = fit->GetParameter(1);
278 const Float_t ndf = fit->GetNDF();
279
280
281 *fLog << inf << "Mean number photo-electrons: " << sig[cnt] << endl;
282 *fLog << inf << "Time Resolution area idx: " << aidx << " Results: " << endl;
283 *fLog << inf << "Mean: " << Form("%4.3f",mean)
284 << "+-" << Form("%4.3f",fit->GetParError(1))
285 << " Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2))
286 << " Chisquare: " << Form("%4.3f",fit->GetChisquare()) << " NDF : " << ndf << endl;
287
288 delete h;
289 gROOT->GetListOfFunctions()->Remove(fit);
290
291 if (sigma < reserr[cnt] && ndf > 2)
292 {
293 res [cnt] = mean;
294 reserr[cnt] = sigma;
295 }
296
297 res[cnt] *= norm;
298 reserr[cnt] *= norm;
299 }
300 else
301 {
302 res[cnt] = -1.;
303 reserr[cnt] = 0.;
304 }
305 cnt++;
306 }
307
308 TGraphErrors *gr = new TGraphErrors(size,
309 sig.GetArray(),res.GetArray(),
310 sigerr.GetArray(),reserr.GetArray());
311 gr->SetTitle(Form("%s%3i","Area Index ",aidx));
312 gr->GetXaxis()->SetTitle("<phes> [1]");
313 gr->GetYaxis()->SetTitle("Time Resolution [ns]");
314 return gr;
315}
Note: See TracBrowser for help on using the repository browser.