source: branches/MarsMoreSimulationTruth/mcalib/MCalibrationIntensityRelTimeCam.cc

Last change on this file was 6680, checked in by gaug, 20 years ago
*** empty log message ***
File size: 15.2 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, const MCalibrationCam::PulserColor_t col)
86{
87
88 if (chargecam.GetSize() != GetSize())
89 {
90 *fLog << err << GetDescriptor() << ": Size mismatch between MCalibrationIntensityRelTimeCam "
91 << "and MCalibrationIntensityChargeCam. " << endl;
92 return NULL;
93 }
94
95 Int_t size = CountNumEntries(col);
96
97 if (size == 0)
98 return NULL;
99
100 if (size != chargecam.CountNumEntries(col))
101 {
102 *fLog << err << GetDescriptor() << ": Size mismatch in colour between MCalibrationIntensityRelTimeCam "
103 << "and MCalibrationIntensityChargeCam. " << endl;
104 return NULL;
105 }
106
107 const Int_t nvalid = chargecam.CountNumValidEntries(pixid,col);
108
109 if (nvalid == 0)
110 {
111 *fLog << err << GetDescriptor() << ": Only un-calibrated events in pixel: " << pixid << endl;
112 return NULL;
113 }
114
115 TArrayF res(nvalid);
116 TArrayF reserr(nvalid);
117 TArrayF sig(nvalid);
118 TArrayF sigerr(nvalid);
119
120 const Float_t sqrt2 = 1.414;
121 const Float_t fadc2ns = 3.333;
122
123 Int_t cnt = 0;
124
125 for (Int_t i=0;i<GetSize();i++)
126 {
127 //
128 // Get the calibration cam from the intensity cam
129 //
130 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)chargecam.GetCam(i);
131 MCalibrationRelTimeCam *relcam = (MCalibrationRelTimeCam*)GetCam(i);
132
133 if (col != MCalibrationCam::kNONE)
134 if (relcam->GetPulserColor() != col)
135 continue;
136 //
137 // Get the calibration pix from the calibration cam
138 //
139 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[pixid];
140 MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[pixid];
141 //
142 // Don't use bad pixels
143 //
144 if (!pix.IsFFactorMethodValid())
145 continue;
146 //
147 res[cnt] = relpix.GetTimePrecision() / sqrt2 * fadc2ns;
148 reserr[cnt] = relpix.GetTimePrecisionErr() / sqrt2 * fadc2ns;
149 //
150 sig [cnt] = pix.GetPheFFactorMethod();
151 sigerr[cnt] = pix.GetPheFFactorMethodErr();
152 cnt++;
153 }
154
155 TGraphErrors *gr = new TGraphErrors(nvalid,
156 sig.GetArray(),res.GetArray(),
157 sigerr.GetArray(),reserr.GetArray());
158 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
159 gr->GetXaxis()->SetTitle("<Photo-electrons> [1]");
160 gr->GetYaxis()->SetTitle("Time Resolution [ns]");
161 return gr;
162}
163
164
165TGraphErrors *MCalibrationIntensityRelTimeCam::GetTimeResVsChargePerArea( const Int_t aidx,const MCalibrationIntensityChargeCam &chargecam, const MGeomCam &geom, const MCalibrationCam::PulserColor_t col)
166{
167
168 Int_t size = CountNumEntries(col);
169
170 if (size == 0)
171 return NULL;
172
173 if (size != chargecam.CountNumEntries(col))
174 {
175 *fLog << err << GetDescriptor() << ": Size mismatch in colour between MCalibrationIntensityRelTimeCam "
176 << "and MCalibrationIntensityChargeCam. " << endl;
177 return NULL;
178 }
179
180 if (col == MCalibrationCam::kBLUE)
181 size -= 5;
182 if (col == MCalibrationCam::kNONE)
183 size -= 5;
184
185 const Float_t sqrt2 = 1.414;
186 const Float_t fadc2ns = 3.333;
187 const Float_t norm = fadc2ns / sqrt2;
188
189 TArrayD res(size);
190 TArrayD reserr(size);
191 TArrayD sig(size);
192 TArrayD sigerr(size);
193
194 Int_t cnt = 0;
195
196 TH1D *h = 0;
197
198 for (Int_t i=0;i<GetSize();i++)
199 {
200 //
201 // Get the calibration cam from the intensity cam
202 //
203 MCalibrationRelTimeCam *relcam = (MCalibrationRelTimeCam*)GetCam(i);
204 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)chargecam.GetCam(i);
205
206 if (relcam->GetPulserColor() != col && col != MCalibrationCam::kNONE)
207 continue;
208
209 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
210 const Double_t phe = (Double_t)apix.GetPheFFactorMethod();
211 const Double_t pheerr = (Double_t)apix.GetPheFFactorMethodErr();
212
213 if (relcam->GetPulserColor() == MCalibrationCam::kBLUE)
214 if (aidx == 0)
215 {
216 if (phe > 100. && phe < 190.)
217 continue;
218 }
219 else
220 {
221 if (phe > 200. && phe < 480.)
222 continue;
223 }
224
225 if (relcam->GetPulserColor() == MCalibrationCam::kUV)
226 *fLog << inf << "NUMBER: " << i << endl;
227
228 sig[cnt] = phe;
229 sigerr[cnt] = pheerr;
230
231 Double_t resol = 0.;
232 Double_t resol2 = 0.;
233 Double_t var = 0.;
234 Int_t num = 0;
235
236 MHCamera camres(geom,"CamRes","Time Resolution;Time Resolution [ns];channels");
237 //
238 // Get the area calibration pix from the calibration cam
239 //
240 for (Int_t j=0; j<cam->GetSize(); j++)
241 {
242 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[j];
243 const MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[j];
244 //
245 // Don't use bad pixels
246 //
247 if (!pix.IsFFactorMethodValid())
248 continue;
249 //
250 //
251 if (aidx != geom[j].GetAidx())
252 continue;
253
254 const Double_t pres = (Double_t)relpix.GetTimePrecision();
255
256 resol += pres;
257 resol2 += pres;
258 num++;
259
260 camres.Fill(j,pres);
261 camres.SetUsed(j);
262 }
263
264 if (num > 100)
265 {
266 resol /= num;
267 var = (resol2 - resol*resol*num) / (num-1);
268
269 res[cnt] = resol;
270 if (var > 0.)
271 reserr[cnt] = TMath::Sqrt(var);
272 else
273 reserr[cnt] = 0.;
274
275 //
276 // Make also a Gauss-fit to the distributions. The RMS can be determined by
277 // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
278 //
279 h = camres.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
280 h->SetDirectory(NULL);
281 h->Fit("gaus","QL");
282 TF1 *fit = h->GetFunction("gaus");
283
284 Double_t ci2 = fit->GetChisquare();
285 Double_t sigma = fit->GetParameter(2);
286
287 if (ci2 > 500. || sigma > reserr[cnt])
288 {
289 h->Fit("gaus","QLM");
290 fit = h->GetFunction("gaus");
291
292 ci2 = fit->GetChisquare();
293 sigma = fit->GetParameter(2);
294 }
295
296 const Double_t mean = fit->GetParameter(1);
297 const Int_t ndf = fit->GetNDF();
298
299 *fLog << inf << "Mean number photo-electrons: " << sig[cnt] << endl;
300 *fLog << inf << "Time Resolution area idx: " << aidx << " Results: " << endl;
301 *fLog << inf << "Mean: " << Form("%4.3f",mean)
302 << "+-" << Form("%4.3f",fit->GetParError(1))
303 << " Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2))
304 << " Chisquare: " << Form("%4.3f",fit->GetChisquare()) << " NDF : " << ndf << endl;
305
306 delete h;
307 gROOT->GetListOfFunctions()->Remove(fit);
308
309 if ((sigma < reserr[cnt] || reserr[cnt]<0.001) && ndf > 2)
310 {
311 res [cnt] = mean;
312 reserr[cnt] = sigma;
313 }
314 else
315 *fLog << warn << "Do not take fit results, but Mean and RMS: " << Form("%4.3f",res[cnt]) << "+-" << Form("%4.3f",reserr[cnt]) << endl;
316
317 res[cnt] *= norm;
318 reserr[cnt] *= norm;
319 cnt++;
320 }
321 }
322
323 TGraphErrors *gr = new TGraphErrors(size,
324 sig.GetArray(),res.GetArray(),
325 sigerr.GetArray(),reserr.GetArray());
326 gr->SetTitle(Form("%s%3i","Area Index ",aidx));
327 gr->GetXaxis()->SetTitle("<phes> [1]");
328 gr->GetYaxis()->SetTitle("Time Resolution [ns]");
329 return gr;
330}
331
332TGraphErrors *MCalibrationIntensityRelTimeCam::GetTimeResVsSqrtPhePerArea( const Int_t aidx,const MCalibrationIntensityChargeCam &chargecam, const MGeomCam &geom, const MCalibrationCam::PulserColor_t col)
333{
334
335 const Int_t size = CountNumEntries(col);
336
337 if (size == 0)
338 return NULL;
339
340 const Int_t asiz = chargecam.CountNumEntries(col);
341
342 if (size != asiz)
343 {
344 *fLog << err << GetDescriptor() << ": Size mismatch in colour between MCalibrationIntensityRelTimeCam "
345 << "and MCalibrationIntensityChargeCam. " << endl;
346 return NULL;
347 }
348
349 const Float_t sqrt2 = 1.414;
350 const Float_t fadc2ns = 3.333;
351 const Float_t norm = fadc2ns / sqrt2;
352
353 TArrayF res(size);
354 TArrayF reserr(size);
355 TArrayF sig(size);
356 TArrayF sigerr(size);
357
358 Int_t cnt = 0;
359
360 TH1D *h = 0;
361
362 for (Int_t i=0;i<GetSize();i++)
363 {
364 //
365 // Get the calibration cam from the intensity cam
366 //
367 MCalibrationRelTimeCam *relcam = (MCalibrationRelTimeCam*)GetCam(i);
368 MCalibrationChargeCam *cam = (MCalibrationChargeCam*)chargecam.GetCam(i);
369
370 if (relcam->GetPulserColor() != col && col != MCalibrationCam::kNONE)
371 continue;
372
373 const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
374 const Float_t phe = apix.GetPheFFactorMethod();
375 const Float_t phesq = phe > 0. ? TMath::Sqrt(phe) : 0.;
376 const Float_t phesqerr = phe > 0. ? 0.5*apix.GetPheFFactorMethodErr()/phesq : 0.;
377
378 sig[cnt] = phesq;
379 sigerr[cnt] = phesqerr;
380
381 Double_t resol = 0.;
382 Double_t resol2 = 0.;
383 Double_t var = 0.;
384 Int_t num = 0;
385
386 MHCamera camres(geom,"CamRes","Time Resolution;Time Resolution [ns];channels");
387 //
388 // Get the area calibration pix from the calibration cam
389 //
390 for (Int_t j=0; j<cam->GetSize(); j++)
391 {
392 const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[j];
393 const MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[j];
394 //
395 // Don't use bad pixels
396 //
397 if (!pix.IsFFactorMethodValid())
398 continue;
399 //
400 //
401 if (aidx != geom[j].GetAidx())
402 continue;
403
404 const Float_t pres = relpix.GetTimePrecision();
405
406 resol += pres;
407 resol2 += pres;
408 num++;
409
410 camres.Fill(j,pres);
411 camres.SetUsed(j);
412 }
413
414 if (num > 1)
415 {
416 resol /= num;
417 var = (resol2 - resol*resol*num) / (num-1);
418
419 res[cnt] = resol;
420 if (var > 0.)
421 reserr[cnt] = TMath::Sqrt(var);
422 else
423 reserr[cnt] = 0.;
424
425 //
426 // Make also a Gauss-fit to the distributions. The RMS can be determined by
427 // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
428 //
429 h = camres.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
430 h->SetDirectory(NULL);
431 h->Fit("gaus","QL");
432 TF1 *fit = h->GetFunction("gaus");
433
434 Float_t ci2 = fit->GetChisquare();
435 Float_t sigma = fit->GetParameter(2);
436
437 if (ci2 > 500. || sigma > reserr[cnt])
438 {
439 h->Fit("gaus","QLM");
440 fit = h->GetFunction("gaus");
441
442 ci2 = fit->GetChisquare();
443 sigma = fit->GetParameter(2);
444 }
445
446 const Float_t mean = fit->GetParameter(1);
447 const Float_t ndf = fit->GetNDF();
448
449
450 *fLog << inf << "Sqrt Mean number photo-electrons: " << sig[cnt] << endl;
451 *fLog << inf << "Time Resolution area idx: " << aidx << " Results: " << endl;
452 *fLog << inf << "Mean: " << Form("%4.3f",mean)
453 << "+-" << Form("%4.3f",fit->GetParError(1))
454 << " Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2))
455 << " Chisquare: " << Form("%4.3f",fit->GetChisquare()) << " NDF : " << ndf << endl;
456
457 delete h;
458 gROOT->GetListOfFunctions()->Remove(fit);
459
460 if (sigma < reserr[cnt] && ndf > 2)
461 {
462 res [cnt] = mean;
463 reserr[cnt] = sigma;
464 }
465
466 res[cnt] *= norm;
467 reserr[cnt] *= norm;
468 }
469 else
470 {
471 res[cnt] = -1.;
472 reserr[cnt] = 0.;
473 }
474 cnt++;
475 }
476
477 TGraphErrors *gr = new TGraphErrors(size,
478 sig.GetArray(),res.GetArray(),
479 sigerr.GetArray(),reserr.GetArray());
480 gr->SetTitle(Form("%s%3i","Area Index ",aidx));
481 gr->GetXaxis()->SetTitle("Sqrt(<phes>) [1]");
482 gr->GetYaxis()->SetTitle("Time Resolution [ns]");
483 return gr;
484}
Note: See TracBrowser for help on using the repository browser.