source: trunk/MagicSoft/Mars/mcalib/MCalibCalcFromPast.cc@ 8199

Last change on this file since 8199 was 8147, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 15.8 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, 12/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24//////////////////////////////////////////////////////////////////////////////
25//
26// MCalibCalcFromPast
27//
28// Steers the occurrance of interlaced calibration events in one data run
29//
30// Input Containers:
31// MParList
32// MCalibrationIntensityChargeCam
33// MCalibrationIntensityRelTimeCam
34// MCalibrationIntensityConstCam
35// MBadPixelsIntensityCam
36//
37// Output Containers:
38// MCalibrationIntensityChargeCam
39// MCalibrationIntensityRelTimeCam
40// MBadPixelsIntensityCam
41//
42// Class version 2:
43// + UInt_t fNumPhesDump; // Number of cams after which the number of phes gets averaged
44// + Float_t fMeanPhes;
45// + Float_t fMeanPhesRelVar;
46// + Bool_t fUpdateNumPhes; // Update the number of photo-electrons only after fNumPhesDump number of Cams
47// + TArrayF fPhes;
48// + TArrayF fPhesVar;
49//
50//////////////////////////////////////////////////////////////////////////////
51#include "MCalibCalcFromPast.h"
52
53#include "MLog.h"
54#include "MLogManip.h"
55
56#include "MParList.h"
57
58#include "MRawRunHeader.h"
59
60#include "MHCalibrationCam.h"
61
62#include "MCalibrationIntensityChargeCam.h"
63#include "MCalibrationIntensityBlindCam.h"
64#include "MCalibrationIntensityQECam.h"
65#include "MCalibrationIntensityRelTimeCam.h"
66#include "MCalibrationIntensityConstCam.h"
67
68#include "MBadPixelsIntensityCam.h"
69
70#include "MCalibrationChargePix.h"
71#include "MCalibrationChargeCalc.h"
72#include "MCalibrationRelTimeCalc.h"
73#include "MCalibrateData.h"
74
75ClassImp(MCalibCalcFromPast);
76
77using namespace std;
78
79const UInt_t MCalibCalcFromPast::fgNumEventsDump = 500;
80const UInt_t MCalibCalcFromPast::fgNumPhesDump = 10;
81
82// --------------------------------------------------------------------------
83//
84// Default constructor.
85//
86// Sets:
87// - fNumEventsDump to fgNumEventsDump
88// - fNumPhesDump to fgNumPhesDump
89//
90MCalibCalcFromPast::MCalibCalcFromPast(const char *name, const char *title)
91 : fGeom(NULL), fParList(NULL), fRunHeader(NULL),
92 fIntensCharge(NULL), fIntensBlind(NULL), fIntensRelTime(NULL), fIntensConst(NULL),
93 fIntensBad(NULL),
94 fChargeCalc(NULL), fRelTimeCalc(NULL), fCalibrate(NULL),
95 fNumCam(0), fNumEvents(0), fUpdateWithFFactorMethod(kTRUE), fUpdateNumPhes(kTRUE),
96 fNumFails(0)
97{
98
99 fName = name ? name : "MCalibCalcFromPast";
100 fTitle = title ? title : "Task to steer the processing of interlaced calibration events";
101
102 SetNumEventsDump();
103 SetNumPhesDump ();
104}
105
106// -----------------------------------------------------------------------------------
107//
108Int_t MCalibCalcFromPast::PreProcess(MParList *pList)
109{
110
111 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
112 if (!fGeom)
113 {
114 *fLog << err << "MGeomCam not found... abort." << endl;
115 return kFALSE;
116 }
117
118 fParList = pList;
119 if (!fParList)
120 {
121 *fLog << err << "MParList not found... abort." << endl;
122 return kFALSE;
123 }
124
125 fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
126 if (!fRunHeader)
127 {
128 *fLog << err << "MRawRunHeader not found... abort." << endl;
129 return kFALSE;
130 }
131
132 //
133 // Look for the MBadPixels Intensity Cam
134 //
135 fIntensBad = (MBadPixelsIntensityCam*)pList->FindCreateObj("MBadPixelsIntensityCam");
136 if (fIntensBad)
137 *fLog << inf << "Found MBadPixelsIntensityCam ... " << flush;
138 else
139 return kFALSE;
140
141 //
142 // Look for the MFillH-class "MHCalibrationBlindCam". In case yes, initialize the
143 // corresponding IntensityCam
144 //
145 if (pList->FindObject(AddSerialNumber("MHCalibrationBlindCam")))
146 {
147
148 *fLog << inf << "Found MHCalibrationBlindCam ... " << flush;
149
150 fIntensBlind = (MCalibrationIntensityBlindCam*)pList->FindCreateObj("MCalibrationIntensityBlindCam");
151 if (!fIntensBlind)
152 {
153 *fLog << err << "Could not find nor create MCalibrationIntensitBlindCam abort... " << endl;
154 return kFALSE;
155 }
156 }
157
158 //
159 // Look for the MFillH-class "MHCalibrationChargeCam". In case yes, initialize the
160 // corresponding IntensityCam
161 //
162 if (pList->FindObject(AddSerialNumber("MHCalibrationChargeCam")))
163 {
164
165 fIntensCharge = (MCalibrationIntensityChargeCam*)pList->FindCreateObj("MCalibrationIntensityChargeCam");
166 fIntensQE = (MCalibrationIntensityQECam*) pList->FindCreateObj("MCalibrationIntensityQECam");
167 fIntensConst = (MCalibrationIntensityConstCam*) pList->FindCreateObj("MCalibrationIntensityConstCam");
168
169 if (!fIntensCharge)
170 return kFALSE;
171 if (!fIntensQE)
172 return kFALSE;
173 if (!fIntensConst)
174 return kFALSE;
175
176
177 MCalibrationChargeCam *chargeinit = (MCalibrationChargeCam*)pList->FindObject("MCalibrationChargeCam");
178 MCalibrationQECam *qeinit = (MCalibrationQECam*) pList->FindObject("MCalibrationQECam");
179
180 if (chargeinit)
181 fIntensCharge->SetCam(chargeinit,0);
182 else
183 *fLog << "Could not find initial MCalibrationChargeCam, cannot initialize intensity cam" << endl;
184
185 if (qeinit)
186 fIntensQE->SetCam(qeinit,0);
187 else
188 *fLog << "Could not find initial MCalibrationQECam, cannot initialize intensity cam" << endl;
189
190 fIntensConst->GetCam()->Init(*fGeom);
191
192 if (!fChargeCalc)
193 fChargeCalc = (MCalibrationChargeCalc*)pList->FindObject("MCalibrationChargeCalc");
194
195 if (!fCalibrate)
196 fCalibrate = (MCalibrateData*)pList->FindObject("MCalibrateData");
197
198 *fLog << inf << "Found MHCalibrationChargeCam ... " << flush;
199
200 if (!fChargeCalc)
201 {
202 *fLog << err << "Could not find MCalibrationChargeCalc abort... " << endl;
203 return kFALSE;
204 }
205
206 if (!fCalibrate)
207 {
208 *fLog << err << "Could not find MCalibrateData abort... " << endl;
209 return kFALSE;
210 }
211 }
212
213 //
214 // Look for the MFillH name "FillRelTimeCam". In case yes, initialize the
215 // corresponding IntensityCam
216 //
217 if (pList->FindObject(AddSerialNumber("MHCalibrationRelTimeCam")))
218 {
219
220 fIntensRelTime = (MCalibrationIntensityRelTimeCam*)pList->FindCreateObj("MCalibrationIntensityRelTimeCam");
221 if (!fRelTimeCalc)
222 fRelTimeCalc = (MCalibrationRelTimeCalc*)pList->FindObject(AddSerialNumber("MCalibrationRelTimeCalc"));
223
224 *fLog << inf << "Found MHCalibrationRelTimeCam ... " << flush;
225
226 if (!fIntensRelTime)
227 {
228 *fLog << err << "Could not find nor create MCalibrationIntensityRelTimeCam abort... " << endl;
229 return kFALSE;
230 }
231
232 if (!fRelTimeCalc)
233 {
234 *fLog << err << "Could not find MCalibrationRelTimeCalc abort... " << endl;
235 return kFALSE;
236 }
237 }
238
239 fNumCam = 0;
240 fNumEvents = 0;
241 fNumPhes = 0;
242
243 fChargeCalc->SetUseExternalNumPhes(kFALSE);
244
245 if (fUpdateNumPhes)
246 {
247 fPhes.Set(fNumPhesDump);
248 fPhesVar.Set(fNumPhesDump);
249 }
250
251 fNumFails = 0;
252
253 return kTRUE;
254}
255
256// --------------------------------------------------------------------------
257//
258// Set fNumEvents=0
259//
260// This is necessary because the calibration histograms do reset themself
261// if ReInit is called, so they are empty. MCalibCalcFromPast errornously
262// ignores how many events are in the histograms but count the number
263// itself.
264//
265Bool_t MCalibCalcFromPast::ReInit(MParList *pList)
266{
267 if (fNumEvents>0)
268 *fLog << inf << fNumEvents << " calibration events at the end of the last file have been skipped." << endl;
269
270 fNumEvents = 0;
271
272 return kTRUE;
273}
274
275// --------------------------------------------------------------------------
276//
277// - Initializes new containers in the
278// - Intensity Cams, if the number of calibration events has reach fNumEventsDump.
279// - Executes Finalize() of the MCalibration*Calc classes in that case.
280// - Sets the latest MCalibrationChargeCam as update class into MCalibrateData
281// - Initialize new MCalibration*Cams into the intensity cams.
282//
283Int_t MCalibCalcFromPast::Process()
284{
285 if (fNumEvents++ < fNumEventsDump)
286 return kTRUE;
287
288 fNumEvents = 0;
289 ReInitialize();
290
291 //
292 // Finalize Possible calibration histogram classes...
293 //
294 *fLog << inf << GetDescriptor() << ": Finalize calibration histograms: " << flush;
295
296 // This fills the IntensityCam which are newly created by
297 // ReInitialize with the result of the last calib cycle
298 if (Finalize("MHCalibrationChargeCam")) *fLog << "MHCalibrationChargeCam..." << flush;
299 if (Finalize("MHCalibrationChargeBlindCam")) *fLog << "MHCalibrationChargeBlindCam..." << flush;
300 if (Finalize("MHCalibrationRelTimeCam")) *fLog << "MHCalibrationRelTimeCam..." << flush;
301
302 //
303 // Finalize possible calibration calculation tasks
304 //
305 *fLog << endl;
306 *fLog << inf << "Finalize calibration calculations..." << endl;
307 if (fChargeCalc)
308 {
309 // Finalized Pedestals, Charges, Bad Pixels and all QE cams
310 if (!fChargeCalc->Finalize())
311 {
312 fNumFails++;
313 *fLog << warn << "WARNING - Finalization of charges failed the " << fNumFails << ". time..." << endl;
314 return kTRUE;
315 }
316
317 if (fUpdateNumPhes)
318 {
319 MCalibrationChargePix &avpix =(MCalibrationChargePix&)fIntensCharge->GetCam()->GetAverageArea(0);
320 fPhes [fNumPhes] = avpix.GetPheFFactorMethod();
321 fPhesVar[fNumPhes] = avpix.GetPheFFactorMethodVar();
322
323 fNumPhes++;
324
325 if (fNumPhes == fNumPhesDump)
326 {
327 fNumPhes = 0;
328 if (!UpdateMeanPhes())
329 {
330 *fLog << warn << "Could not update mean number of photo-electrons. "
331 << "Skip it until next update" << endl;
332 fChargeCalc->SetUseExternalNumPhes(kFALSE);
333 }
334 else
335 {
336 *fLog << inf << "New averaged number photo-electrons: " << fMeanPhes << endl;
337 fChargeCalc->SetExternalNumPhes ( fMeanPhes );
338 fChargeCalc->SetExternalNumPhesRelVar( fMeanPhesRelVar );
339 fChargeCalc->SetUseExternalNumPhes();
340 }
341 }
342 }
343 }
344
345 if (fRelTimeCalc)
346 fRelTimeCalc->Finalize();
347
348 if (fCalibrate)
349 fCalibrate->UpdateConversionFactors(fUpdateWithFFactorMethod ? NULL
350 : (MCalibrationChargeCam*)fIntensCharge->GetCam() );
351
352 return kTRUE;
353}
354
355
356// --------------------------------------------------------------------------
357//
358// Searches for name in the MParList and calls, if existing:
359// - MHCalibrationCam::Finalize()
360// - MHCalibrationCam::ResetHists()
361//
362Bool_t MCalibCalcFromPast::Finalize(const char* name, Bool_t finalize)
363{
364 MHCalibrationCam *hist = (MHCalibrationCam*)fParList->FindObject(name, "MHCalibrationCam");
365 if (!hist)
366 return kFALSE;
367
368 if (finalize)
369 hist->Finalize();
370
371 hist->ResetHists();
372 return kTRUE;
373}
374
375// --------------------------------------------------------------------------
376//
377// Re-Intitializes new containers inside the Intensity Cams.
378// From now on, a call to the IntensityCam functions returns pointers
379// to the newly created Containers.
380//
381Bool_t MCalibCalcFromPast::ReInitialize()
382{
383 fNumCam++;
384
385 *fLog << inf << "MCalibCalcFromPast::ReInitialize #" << fNumCam << ": ";
386
387 const Int_t runnumber = fRunHeader->GetRunNumber();
388
389 // The "DeleteOldCam" function must not delete the first entry in
390 // the array because it is a special cam from the MParList. (see above)
391
392 if (fIntensBad)
393 {
394 *fLog << "MBadPixelsCam...";
395 MBadPixelsCam *cold = fIntensBad->GetCam();
396 fIntensBad->AddToList(Form("MBadPixelsCam%04d",fNumCam),*fGeom);
397 if (cold)
398 fIntensBad->DeleteOldCam(cold);
399 }
400
401 if (fIntensCharge)
402 {
403 *fLog << "MCalibrationChargeCam...";
404 MCalibrationChargeCam *cold = (MCalibrationChargeCam*)fIntensCharge->GetCam();
405
406 fIntensCharge->AddToList(Form("MCalibrationChargeCam%04d",fNumCam),*fGeom);
407
408 MCalibrationChargeCam *cnew = (MCalibrationChargeCam*)fIntensCharge->GetCam();
409 cnew->SetRunNumber(runnumber);
410
411 if (cold)
412 {
413 cnew->MergeHiLoConversionFactors(*cold);
414 fIntensCharge->DeleteOldCam(cold);
415 }
416 }
417
418 if (fIntensQE)
419 {
420 *fLog << "MCalibrationQECam...";
421 MCalibrationCam *cold = fIntensQE->GetCam();
422 fIntensQE->AddToList(Form("MCalibrationQECam%04d",fNumCam),*fGeom);
423 fIntensQE->GetCam()->SetRunNumber(runnumber);
424 if (cold)
425 fIntensQE->DeleteOldCam(cold);
426 }
427 if (fIntensBlind)
428 {
429 *fLog << "MCalibrationBlindCam...";
430 MCalibrationCam *cold = fIntensBlind->GetCam();
431 fIntensBlind->AddToList(Form("MCalibrationBlindCam%04d",fNumCam),*fGeom);
432 fIntensBlind->GetCam()->SetRunNumber(runnumber);
433 if (cold)
434 fIntensBlind->DeleteOldCam(cold);
435 }
436
437 *fLog << endl;
438
439 return kTRUE;
440
441}
442
443Int_t MCalibCalcFromPast::PostProcess()
444{
445 if (GetNumExecutions()==0)
446 return kTRUE;
447
448 // Now we reset all histograms to make sure that the PostProcess
449 // of the following tasks doesn't try to finalize a partly empty
450 // histogram!
451 Finalize("MHCalibrationChargeCam", kFALSE);
452 Finalize("MHCalibrationChargeBlindCam", kFALSE);
453 Finalize("MHCalibrationRelTimeCam", kFALSE);
454
455 if (fChargeCalc)
456 fChargeCalc->ResetNumProcessed();
457
458 if (fNumCam==0)
459 return kTRUE;
460
461 *fLog << inf << endl;
462 *fLog << GetDescriptor() << " execution statistics:" << endl;
463 *fLog << " " << setfill(' ') << setw(7) << fNumFails << " (" << Form("%5.1f", 100.*fNumFails/fNumCam) << "%) updates failed." << endl;
464 *fLog << endl;
465
466 return kTRUE;
467}
468
469
470Bool_t MCalibCalcFromPast::UpdateMeanPhes()
471{
472 Float_t sumw = 0.;
473 Float_t sum = 0.;
474
475 for (Int_t i=0; i<fPhes.GetSize(); i++)
476 {
477 const Float_t weight = 1./fPhesVar[i];
478 sum += fPhes[i]*weight;
479 sumw += weight;
480 }
481
482 if (sumw < 0.000001)
483 return kFALSE;
484
485 if (sum < 0.000001)
486 return kFALSE;
487
488 fMeanPhes = sum/sumw;
489 fMeanPhesRelVar = sumw/sum/sum;
490
491 return kTRUE;
492}
493
494// --------------------------------------------------------------------------
495//
496// Read the setup from a TEnv, eg:
497// MCalibCalcFromPast.UpdateWithFFactorMethod: Off, On
498// MCalibCalcFromPast.NumEventsDump: 500
499// MCalibCalcFromPast.UpdateNumPhes: yes/no
500// MCalibCalcFromPast.NumPhesDump: 10
501//
502Int_t MCalibCalcFromPast::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
503{
504 Bool_t rc = kFALSE;
505 if (IsEnvDefined(env, prefix, "UpdateWithFFactorMethod", print))
506 {
507 rc = kTRUE;
508 SetUpdateWithFFactorMethod(GetEnvValue(env, prefix, "UpdateWithFFactorMethod", fUpdateWithFFactorMethod));
509 }
510
511 if (IsEnvDefined(env, prefix, "NumEventsDump", print))
512 {
513 SetNumEventsDump(GetEnvValue(env, prefix, "NumEventsDump", (Int_t)fNumEventsDump));
514 rc = kTRUE;
515 }
516
517 if (IsEnvDefined(env, prefix, "UpdateNumPhes", print))
518 {
519 TString s = GetEnvValue(env, prefix, "UpdateNumPhes", "");
520 s.ToLower();
521 if (s.BeginsWith("no"))
522 SetUpdateNumPhes(kFALSE);
523 if (s.BeginsWith("yes"))
524 SetUpdateNumPhes(kTRUE);
525 rc = kTRUE;
526 }
527
528 if (IsEnvDefined(env, prefix, "NumPhesDump", print))
529 {
530 SetNumPhesDump(GetEnvValue(env, prefix, "NumPhesDump", (Int_t)fNumPhesDump));
531 rc = kTRUE;
532 }
533
534 return rc;
535}
Note: See TracBrowser for help on using the repository browser.