source: trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc@ 2643

Last change on this file since 2643 was 2642, checked in by gaug, 22 years ago
*** empty log message ***
File size: 13.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!
19! Author(s): Markus Gaug 11/2003 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2001
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27// //
28// MCalibrationCam //
29// //
30// Hold the Calibration information for all pixels in the camera //
31// //
32// The Classes MCalibrationPix's are stored in a TClonesArray //
33// The Class MCalibrationBlindPix and the MCalibrationPINDiode //
34// are stored in separate pointers //
35// //
36/////////////////////////////////////////////////////////////////////////////
37#include "MCalibrationCam.h"
38#include "MCalibrationPix.h"
39#include "MCalibrationBlindPix.h"
40#include "MCalibrationConfig.h"
41
42#include <TClonesArray.h>
43
44#include "MLog.h"
45#include "MLogManip.h"
46
47#include "MGeomCam.h"
48
49ClassImp(MCalibrationCam);
50
51using namespace std;
52// --------------------------------------------------------------------------
53//
54// Default constructor.
55// Creates a TClonesArray of MHCalibrationPix objects, initialized to 0 entries
56// Creates an MCalibrationBlindPix and an MCalibrationPINDiode
57//
58MCalibrationCam::MCalibrationCam(const char *name, const char *title)
59 : fMeanNrPhotAvailable(kFALSE),
60 fMeanNrPhotInnerPix(-1.),
61 fMeanNrPhotInnerPixErr(-1.)
62{
63 fName = name ? name : "MCalibrationCam";
64 fTitle = title ? title : "Storage container for the Calibration Information in the camera";
65
66 fPixels = new TClonesArray("MCalibrationPix",1);
67 fBlindPixel = new MCalibrationBlindPix();
68 fPINDiode = new MCalibrationPINDiode();
69}
70
71// --------------------------------------------------------------------------
72//
73// Delete the TClonesArray or MCalibrationPix's
74// Delete the MCalibrationPINDiode and the MCalibrationBlindPix
75//
76MCalibrationCam::~MCalibrationCam()
77{
78
79 //
80 // delete fPixels should delete all Objects stored inside
81 //
82 delete fPixels;
83 delete fBlindPixel;
84 delete fPINDiode;
85}
86
87// --------------------------------------------------------------------------
88//
89// This function return the size of the FILLED MCalibrationCam
90// It is NOT the size of the array fPixels !!!
91// Note that with every call to AddPixels, GetSize() might change
92//
93Int_t MCalibrationCam::GetSize() const
94{
95
96 //
97 // Here we get the number of "filled" fPixels!!
98 //
99 return fPixels->GetEntriesFast();
100}
101
102// --------------------------------------------------------------------------
103//
104// Get i-th pixel (pixel number)
105//
106MCalibrationPix &MCalibrationCam::operator[](Int_t i)
107{
108 return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
109}
110
111// --------------------------------------------------------------------------
112//
113// Get i-th pixel (pixel number)
114//
115MCalibrationPix &MCalibrationCam::operator[](Int_t i) const
116{
117 return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
118}
119
120// --------------------------------------------------------------------------
121//
122// Return a pointer to the pixel with the requested idx.
123// NULL if it doesn't exist.
124//
125MCalibrationPix *MCalibrationCam::GetCalibrationPix(Int_t idx) const
126{
127 if (idx<0)
128 return NULL;
129
130 if (!CheckBounds(idx))
131 return NULL;
132
133 return (MCalibrationPix*)fPixels->At(idx);
134}
135
136
137
138// --------------------------------------------------------------------------
139//
140// Check if position i is inside bounds
141//
142Bool_t MCalibrationCam::CheckBounds(Int_t i) const
143{
144 return i < fPixels->GetEntriesFast();
145}
146
147Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const
148{
149 if (!CheckBounds(idx))
150 return kFALSE;
151
152 return kTRUE;
153}
154
155Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const
156{
157 return ((*this)[idx].GetCharge() > 0. && (*this)[idx].GetErrCharge() > 0.);
158}
159
160
161
162void MCalibrationCam::Clear(Option_t *o)
163{
164 fPixels->ForEach(TObject, Clear)();
165}
166
167//
168// Perform the fits of the charges
169// If i=-1, then all pixels will be fitted
170// Otherwise only the one with index i
171//
172// The number of succesful fits is returned
173//
174UShort_t MCalibrationCam::FitCharge(Int_t i)
175{
176
177 UShort_t nsuccess = 0;
178
179 // invalidate i if it exceeds the number of entries in fPixels
180 if (i > fPixels->GetEntriesFast())
181 {
182 *fLog << warn << "Tried to fit pixel out of allowed range " << endl;
183 return 0;
184 }
185
186 if (i == -1) // loop over all events
187 {
188
189 TIter Next(fPixels);
190 MCalibrationPix *pix;
191 while ((pix=(MCalibrationPix*)Next()))
192 {
193 if (pix->FitCharge())
194 nsuccess++;
195 }
196 }
197 else // fit only the pixel with index i
198 {
199 if((*this)[i].FitCharge())
200 nsuccess++;
201 }
202
203 return nsuccess;
204
205}
206
207//
208// Perform the fits of the charges of all pixels
209// plus the blind pixel and the PIN Diode
210//
211// The number of succesful fits is returned
212//
213UShort_t MCalibrationCam::FitAllCharge()
214{
215
216 // FIXME: Once the blind pixel is fully working,
217 // there must be some penalty in case the
218 // fit does not succeed
219 //
220 UShort_t nsuccess = 0;
221
222 TIter Next(fPixels);
223 MCalibrationPix *pix;
224 while ((pix=(MCalibrationPix*)Next()))
225 {
226 if (pix->FitCharge())
227 nsuccess++;
228 }
229
230 if (fBlindPixel->FitCharge())
231 nsuccess++;
232
233 if (fPINDiode->FitCharge())
234 nsuccess++;
235
236 return nsuccess;
237
238}
239
240
241
242//
243// Perform the fits of the arrival times
244// If i=-1, then all pixels will be fitted
245// Otherwise only the one with index i
246//
247// The number of succesful fits is returned
248//
249UShort_t MCalibrationCam::FitTime(Int_t i)
250{
251
252 UShort_t nsuccess = 0;
253
254 // invalidate i if it exceeds the number of entries in fPixels
255 if (i > fPixels->GetEntriesFast())
256 {
257 *fLog << warn << "Tried to fit pixel out of allowed range " << endl;
258 return 0;
259 }
260
261 if (i == -1) // loop over all events
262 {
263
264 TIter Next(fPixels);
265 MCalibrationPix *pix;
266 while ((pix=(MCalibrationPix*)Next()))
267 {
268 if (pix->FitTime())
269 nsuccess++;
270 }
271 }
272 else // fit only the pixel with index i
273 {
274 if((*this)[i].FitTime())
275 nsuccess++;
276 }
277
278 return nsuccess;
279
280}
281
282
283//
284// Perform the fits of the times of all pixels
285// plus the blind pixel and the PIN Diode
286//
287// The number of succesful fits is returned
288//
289UShort_t MCalibrationCam::FitAllTime()
290{
291
292 // FIXME: Once the blind pixel is fully working,
293 // there must be some penalty in case the
294 // fit does not succeed
295 //
296 UShort_t nsuccess = 0;
297
298 TIter Next(fPixels);
299 MCalibrationPix *pix;
300 while ((pix=(MCalibrationPix*)Next()))
301 {
302 if (pix->FitTime())
303 nsuccess++;
304 }
305
306 if (fBlindPixel->FitTime())
307 nsuccess++;
308
309 if (fPINDiode->FitTime())
310 nsuccess++;
311
312 return nsuccess;
313
314}
315
316
317
318void MCalibrationCam::CutEdges()
319{
320
321 fBlindPixel->GetHist()->CutAllEdges();
322
323 TIter Next(fPixels);
324 MCalibrationPix *pix;
325 while ((pix=(MCalibrationPix*)Next()))
326 {
327 pix->GetHist()->CutAllEdges();
328 }
329
330 return;
331}
332
333// ---------------------------------------------------------------------
334//
335// This function simply allocates memory via the ROOT command:
336// (TObject**) TStorage::ReAlloc(fCont, newSize * sizeof(TObject*),
337// fSize * sizeof(TObject*));
338// newSize corresponds to size in our case
339// fSize is the old size (in most cases: 0)
340//
341void MCalibrationCam::InitSize(Int_t size)
342{
343
344 //
345 // check if we have already initialized to size
346 //
347 if (CheckBounds(size))
348 return;
349
350 //
351 // It is important to use Expand and NOT ExpandCreate, because
352 // we want to keep all pixel not used with a NULL pointer.
353 //
354 fPixels->ExpandCreate(size);
355
356}
357
358void MCalibrationCam::Print(Option_t *o) const
359{
360 *fLog << all << GetDescriptor() << ":" << endl;
361 int id = 0;
362
363 *fLog << "Succesfully calibrated pixels:" << endl;
364 *fLog << endl;
365
366 TIter Next(fPixels);
367 MCalibrationPix *pix;
368 while ((pix=(MCalibrationPix*)Next()))
369 {
370
371 if (pix->GetCharge() >= 0.)
372 {
373 *fLog << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms()
374 << " Reduced Charge: " << pix->GetCharge() << " +- "
375 << pix->GetSigmaCharge() << " Reduced Sigma: " << TMath::Sqrt(pix->GetRSigmaSquare()) << endl;
376 id++;
377 }
378 }
379
380 *fLog << id << " succesful pixels :-))" << endl;
381 id = 0;
382
383 *fLog << endl;
384 *fLog << "Pixels with errors:" << endl;
385 *fLog << endl;
386
387 TIter Next2(fPixels);
388 while ((pix=(MCalibrationPix*)Next2()))
389 {
390
391 if (pix->GetCharge() == -1.)
392 {
393 *fLog << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms()
394 << " Reduced Charge: " << pix->GetCharge() << " +- "
395 << pix->GetSigmaCharge() << " Reduced Sigma: " << TMath::Sqrt(pix->GetRSigmaSquare()) << endl;
396 id++;
397 }
398 }
399 *fLog << id << " pixels with errors :-((" << endl;
400
401}
402
403
404Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
405{
406
407 switch (type)
408 {
409 case 0:
410 val = (*this)[idx].GetCharge();
411 break;
412 case 1:
413 val = (*this)[idx].GetErrCharge();
414 break;
415 case 2:
416 val = (*this)[idx].GetSigmaCharge();
417 break;
418 case 3:
419 val = (*this)[idx].GetErrSigmaCharge();
420 break;
421 case 4:
422 val = (*this)[idx].GetChargeProb();
423 break;
424 case 5:
425 val = (*this)[idx].GetTime();
426 break;
427 case 6:
428 val = (*this)[idx].GetSigmaTime();
429 break;
430 case 7:
431 val = (*this)[idx].GetTimeProb();
432 break;
433 case 8:
434 val = (*this)[idx].GetPed();
435 break;
436 case 9:
437 val = (*this)[idx].GetPedRms();
438 break;
439 case 10:
440 val = TMath::Sqrt((*this)[idx].GetRSigmaSquare());
441 break;
442 case 15:
443 val = ((*this)[idx].GetSigmaCharge()/(*this)[idx].GetCharge())*
444 ((*this)[idx].GetSigmaCharge()/(*this)[idx].GetCharge());
445 break;
446 case 11:
447 val = (*this)[idx].GetPheFFactorMethod();
448 break;
449 case 12:
450 val = (*this)[idx].GetConversionFFactorMethod();
451 break;
452 case 13:
453 if (idx < 397)
454 val = (double)fMeanNrPhotInnerPix;
455 else
456 val = (double)fMeanNrPhotInnerPix*gkCalibrationOuterPixelArea;
457 break;
458 case 14:
459 if ((fMeanNrPhotInnerPix > 0. ) && ((*this)[idx].GetCharge() != -1.))
460 {
461 if (idx < 397)
462 val = fMeanNrPhotInnerPix / (*this)[idx].GetCharge();
463 else
464 val = fMeanNrPhotInnerPix*gkCalibrationOuterPixelArea / (*this)[idx].GetCharge();
465 }
466 else
467 {
468 val = -1.;
469 }
470 break;
471 default:
472 return kFALSE;
473 }
474 return val>=0;
475}
476
477void MCalibrationCam::DrawPixelContent(Int_t idx) const
478{
479
480 (*this)[idx].Draw();
481
482}
483
484Bool_t MCalibrationCam::CalcNrPhotInnerPixel()
485{
486 if (!fBlindPixel->IsFitOK())
487 return kFALSE;
488
489 const Float_t mean = fBlindPixel->GetLambda();
490 // const Float_t merr = fBlindPixel->GetErrLambda();
491
492 switch (fColor)
493 {
494 case kECGreen:
495 fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQEGreen) // real photons
496 *TMath::Power(10,gkCalibrationBlindPixelAttGreen) // correct for absorption
497 / gkCalibrationBlindPixelArea; // correct for area
498 break;
499 case kECBlue:
500 fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQEBlue )
501 *TMath::Power(10,gkCalibrationBlindPixelAttBlue)
502 / gkCalibrationBlindPixelArea;
503 break;
504 case kECUV:
505 fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQEUV )
506 *TMath::Power(10,gkCalibrationBlindPixelAttUV)
507 / gkCalibrationBlindPixelArea;
508 break;
509 case kECCT1:
510 default:
511 fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQECT1 )
512 *TMath::Power(10,gkCalibrationBlindPixelAttCT1)
513 / gkCalibrationBlindPixelArea;
514 break;
515 }
516
517 fMeanNrPhotAvailable = kTRUE;
518 return kTRUE;
519}
520
521
522Bool_t MCalibrationCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
523{
524
525 if (ipx < 0 || !IsPixelFitted(ipx))
526 return kFALSE;
527
528 if (!fMeanNrPhotAvailable)
529 if (!CalcNrPhotInnerPixel())
530 return kFALSE;
531
532 mean = fMeanNrPhotInnerPix / (*this)[ipx].GetCharge();
533
534 //
535 // Not yet ready , sorry
536 //
537 err = -1.;
538 sigma = -1.;
539
540 return kTRUE;
541}
542
543
544Bool_t MCalibrationCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
545{
546
547 if (ipx < 0 || !IsPixelFitted(ipx))
548 return kFALSE;
549
550 Float_t conv = (*this)[ipx].GetConversionFFactorMethod();
551
552 if (conv < 0.)
553 return kFALSE;
554
555 mean = conv;
556
557 //
558 // Not yet ready , sorry
559 //
560 err = -1.;
561 sigma = -1.;
562
563 return kTRUE;
564}
565
566
567Bool_t MCalibrationCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
568{
569
570 return kTRUE;
571
572}
573
574Bool_t MCalibrationCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
575{
576
577 return kTRUE;
578
579}
Note: See TracBrowser for help on using the repository browser.