source: trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityConstCam.cc@ 7511

Last change on this file since 7511 was 7195, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 11.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 06/2005 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24/////////////////////////////////////////////////////////////////////////////
25//
26// MCalibrationIntensityConstCam
27//
28// Storage container for intensity charge calibration results.
29//
30// Individual MCalibrationConstCam'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, MCalibConstCam,
36// MCalibrationChargePix, MCalibrationChargeCalc, MCalibrationQECam
37// MHCalibrationPix, MHCalibrationCam
38//
39/////////////////////////////////////////////////////////////////////////////
40#include "MCalibrationIntensityConstCam.h"
41#include "MCalibConstCam.h"
42#include "MCalibConstPix.h"
43
44#include "MGeomCam.h"
45#include "MGeomPix.h"
46
47#include "MBadPixelsCam.h"
48#include "MBadPixelsPix.h"
49
50#include "MLogManip.h"
51#include "MHCamera.h"
52
53#include <TH1D.h>
54#include <TF1.h>
55#include <TOrdCollection.h>
56#include <TGraphErrors.h>
57
58ClassImp(MCalibrationIntensityConstCam);
59
60using namespace std;
61// --------------------------------------------------------------------------
62//
63// Default constructor.
64//
65MCalibrationIntensityConstCam::MCalibrationIntensityConstCam(const char *name, const char *title)
66{
67
68 fName = name ? name : "MCalibrationIntensityConstCam";
69 fTitle = title ? title : "Results of the Intensity Calibration";
70
71 fCams = new TOrdCollection;
72 fCams->SetOwner();
73
74 InitSize(1);
75}
76
77// --------------------------------------------------------------------------
78//
79// Deletes the histograms if they exist
80//
81MCalibrationIntensityConstCam::~MCalibrationIntensityConstCam()
82{
83 if (fCams)
84 delete fCams;
85}
86
87// -------------------------------------------------------------------
88//
89// Add MCalibrationConstCam's in the ranges from - to.
90//
91void MCalibrationIntensityConstCam::Add(const UInt_t from, const UInt_t to)
92{
93 for (UInt_t i=from; i<to; i++)
94 fCams->AddAt(new MCalibConstCam,i);
95}
96
97// --------------------------------------------------------------------------
98//
99// Add a new MCalibConstCam to fCams, give it the name "name" and initialize
100// it with geom.
101//
102void MCalibrationIntensityConstCam::AddToList( const char* name, const MGeomCam &geom)
103{
104 InitSize(GetSize()+1);
105 GetCam()->SetName(name);
106 GetCam()->Init(geom);
107}
108
109// -----------------------------------------------------
110//
111// Calls Clear() for all entries fCams
112//
113void MCalibrationIntensityConstCam::Clear(Option_t *o)
114{
115
116 fCams->ForEach(MCalibConstCam, Clear)();
117 return;
118}
119
120// -----------------------------------------------------
121//
122// Calls Print(o) for all entries fCams
123//
124void MCalibrationIntensityConstCam::Print(Option_t *o) const
125{
126 fCams->ForEach(MCalibConstCam, Print)(o);
127}
128
129
130// -------------------------------------------------------------------
131//
132// Initialize the objects inside the TOrdCollection using the
133// virtual function Add().
134//
135// InitSize can only increase the size, but not shrink.
136//
137// It can be called more than one time. New Containers are
138// added only from the current size to the argument i.
139//
140void MCalibrationIntensityConstCam::InitSize(const UInt_t i)
141{
142
143 const UInt_t save = GetSize();
144
145 if (i==save)
146 return;
147
148 if (i>save)
149 Add(save,i);
150}
151
152// -------------------------------------------------------------------
153//
154// If size is still 0, Intialize a first Cam.
155// Calls Init(geom) for all fCams
156//
157void MCalibrationIntensityConstCam::Init(const MGeomCam &geom)
158{
159 if (GetSize() == 0)
160 InitSize(1);
161
162 fCams->ForEach(MCalibConstCam,Init)(geom);
163}
164
165
166
167// --------------------------------------------------------------------------
168//
169// Returns the current size of the TOrdCollection fCams
170// independently if the MCalibrationCam is filled with values or not.
171//
172const Int_t MCalibrationIntensityConstCam::GetSize() const
173{
174 return fCams->GetSize();
175}
176
177// --------------------------------------------------------------------------
178//
179// Get i-th pixel from current camera
180//
181MCalibConstPix &MCalibrationIntensityConstCam::operator[](UInt_t i)
182{
183 return (*GetCam())[i];
184}
185
186// --------------------------------------------------------------------------
187//
188// Get i-th pixel from current camera
189//
190const MCalibConstPix &MCalibrationIntensityConstCam::operator[](UInt_t i) const
191{
192 return (*GetCam())[i];
193}
194
195// --------------------------------------------------------------------------
196//
197// Get i-th camera
198//
199MCalibConstCam *MCalibrationIntensityConstCam::GetCam(Int_t i)
200{
201 return static_cast<MCalibConstCam*>(i==-1 ? fCams->Last() : fCams->At(i));
202}
203
204// --------------------------------------------------------------------------
205//
206// Get i-th camera
207//
208const MCalibConstCam *MCalibrationIntensityConstCam::GetCam(Int_t i) const
209{
210 return static_cast<MCalibConstCam*>(i==-1 ? fCams->Last() : fCams->At(i));
211}
212
213Bool_t MCalibrationIntensityConstCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
214{
215 if ((*fBadPixels)[idx].IsUnsuitable())
216 return kFALSE;
217
218 return GetCam()->GetPixelContent(val,idx,cam,type);
219}
220
221TGraphErrors *MCalibrationIntensityConstCam::GetConvFactorPerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
222{
223
224 const Int_t size = GetSize();
225
226 if (size == 0)
227 return NULL;
228
229 TArrayF convarea(size);
230 TArrayF convareaerr(size);
231 TArrayF time(size);
232 TArrayF timeerr(size);
233
234 TH1D *h = 0;
235 MCalibConstCam *cam = 0;
236
237 for (Int_t i=0;i<GetSize();i++)
238 {
239 //
240 // Get the calibration cam from the intensity cam
241 //
242 cam = GetCam(i);
243
244 //
245 // Get the calibration pix from the calibration cam
246 //
247 Double_t variab = 0.;
248 Double_t variab2 = 0.;
249 Double_t variance = 0.;
250 Int_t num = 0;
251 Float_t pvar = 0.;
252
253 MHCamera camconv(geom,"CamConv","Variable;;channels");
254 //
255 // Get the area calibration pix from the calibration cam
256 //
257 for (Int_t j=0; j<cam->GetSize(); j++)
258 {
259 const MCalibConstPix &pix = (*cam)[j];
260 //
261 // Use only pixels of a same area index
262 //
263 if (aidx != geom[j].GetAidx())
264 continue;
265
266 pvar = pix.GetCalibConst();
267
268 if (pvar < 0)
269 continue;
270
271 variab += pvar;
272 variab2 += pvar*pvar;
273 num++;
274
275 camconv.Fill(j,pvar);
276 camconv.SetUsed(j);
277 }
278
279 if (num > 1)
280 {
281 variab /= num;
282 variance = (variab2 - variab*variab*num) / (num-1);
283
284 convarea[i] = variab;
285 if (variance > 0.)
286 convareaerr[i] = TMath::Sqrt(variance);
287 else
288 convareaerr[i] = 999999999.;
289
290 //
291 // Make also a Gauss-fit to the distributions. The RMS can be determined by
292 // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
293 //
294 h = camconv.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",500);
295 h->SetDirectory(NULL);
296 h->Fit("gaus","QL");
297 TF1 *fit = h->GetFunction("gaus");
298
299 Float_t ci2 = fit->GetChisquare();
300 Float_t sigma = fit->GetParameter(2);
301
302 if (ci2 > 500. || sigma > convareaerr[i])
303 {
304 h->Fit("gaus","QLM");
305 fit = h->GetFunction("gaus");
306
307 ci2 = fit->GetChisquare();
308 sigma = fit->GetParameter(2);
309 }
310
311 const Float_t mean = fit->GetParameter(1);
312 const Float_t meanerr = fit->GetParError(1);
313 const Float_t ndf = fit->GetNDF();
314
315 *fLog << inf << "Camera Nr: " << i << endl;
316 *fLog << inf << " area idx: " << aidx << " Results: " << endl;
317 *fLog << inf << "Mean: " << Form("%4.3f",mean)
318 << "+-" << Form("%4.3f",meanerr)
319 << " Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2))
320 << " Chisquare: " << Form("%4.3f",ci2) << " NDF : " << ndf << endl;
321 delete h;
322 gROOT->GetListOfFunctions()->Remove(fit);
323
324 if (meanerr < convareaerr[i] && ndf > 2 && ci2/ndf < 10.)
325 {
326 convarea [i] = mean;
327 convareaerr[i] = meanerr;
328 }
329 else
330 {
331 convareaerr[i] /= TMath::Sqrt((Float_t)num);
332 }
333 }
334 else
335 {
336 convarea[i] = -1.;
337 convareaerr[i] = 0.;
338 }
339
340 time[i] = i;
341 timeerr[i] = 0.;
342 }
343
344 TGraphErrors *gr = new TGraphErrors(size,
345 time.GetArray(),convarea.GetArray(),
346 timeerr.GetArray(),convareaerr.GetArray());
347 gr->SetTitle(Form("Conversion Factors Area %3i Average",aidx));
348 gr->GetXaxis()->SetTitle("Camera Nr.");
349 gr->GetYaxis()->SetTitle("<N_{phe}>/<Q> [FADC cnt^{-1}]");
350 return gr;
351}
352
353TGraph *MCalibrationIntensityConstCam::GetConvFactorVsTime ( const Int_t pixid )
354{
355
356 const Int_t size = GetSize();
357
358 if (size == 0)
359 return NULL;
360
361 TArrayF time(size);
362 TArrayF conv(size);
363
364 MCalibConstCam *cam = 0;
365
366 for (Int_t i=0;i<size;i++)
367 {
368 //
369 // Get the calibration cam from the intensity cam
370 //
371 cam = GetCam(i);
372
373 if (cam->GetSize()-1 < pixid)
374 {
375 *fLog << err << "Pixel " << pixid << " out of range " << endl;
376 return NULL;
377 }
378
379 //
380 // Get the calibration pix from the calibration cam
381 //
382 MCalibConstPix &pix = (*cam)[pixid];
383 //
384 time[i] = i;
385 conv[i] = pix.GetCalibConst();
386 //
387 }
388
389
390 TGraph *gr = new TGraph(size,time.GetArray(),conv.GetArray());
391 gr->SetTitle(Form("%s%3i","Pixel ",pixid));
392 gr->GetXaxis()->SetTitle("Camera Nr.");
393 gr->GetYaxis()->SetTitle("<N_{phe}>/<Q> [FADC cnts^{-1}]");
394 return gr;
395
396
397}
398
399void MCalibrationIntensityConstCam::DrawConvFactorPerAreaVsTime( const Int_t aidx )
400{
401
402 TGraphErrors *gr = GetConvFactorPerAreaVsTime(aidx,MGeomCamMagic());
403
404 if (gr)
405 {
406 gr->SetBit(kCanDelete);
407 gr->Draw("A*");
408 }
409}
410
411void MCalibrationIntensityConstCam::DrawConvFactorVsTime( const Int_t idx )
412{
413 TGraph *gr = GetConvFactorVsTime(idx);
414 if (gr)
415 {
416 gr->SetBit(kCanDelete);
417 gr->Draw("A*");
418 }
419}
420
Note: See TracBrowser for help on using the repository browser.