source: trunk/MagicSoft/Mars/mjobs/MJPedestal.cc@ 8790

Last change on this file since 8790 was 8788, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 38.7 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): Thomas Bretz, 1/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Markus Gaug, 4/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2007
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MJPedestal
29//
30// Resource file entries are case sensitive!
31//
32// We require at least fMinEvents (def=50) to be processed by the
33// ExtractPedestal-task. If not an error is returned.
34//
35/////////////////////////////////////////////////////////////////////////////
36#include "MJPedestal.h"
37
38// C/C++ includes
39#include <fstream>
40
41// root classes
42#include <TF1.h>
43#include <TEnv.h>
44#include <TFile.h>
45#include <TLine.h>
46#include <TLatex.h>
47#include <TString.h>
48#include <TCanvas.h>
49#include <TSystem.h>
50#include <TLegend.h>
51#include <TPad.h>
52#include <TEnv.h>
53#include <TH2F.h>
54
55// mars core
56#include "MLog.h"
57#include "MLogManip.h"
58
59#include "MTaskEnv.h"
60#include "MSequence.h"
61#include "MRunIter.h"
62#include "MParList.h"
63#include "MTaskList.h"
64#include "MEvtLoop.h"
65
66#include "MStatusDisplay.h"
67
68// Other basic classes
69#include "MExtractTimeAndCharge.h"
70
71// parameter containers
72#include "MGeomCam.h"
73#include "MHCamera.h"
74#include "MPedestalPix.h"
75
76#include "MCalibrationPix.h"
77#include "MHPedestalPix.h"
78#include "MHCalibrationPulseTimeCam.h"
79#include "MCalibrationPulseTimeCam.h"
80
81// tasks
82#include "MReadMarsFile.h"
83#include "MRawFileRead.h"
84#include "MRawEvtData.h"
85#include "MGeomApply.h"
86#include "MPedestalSubtract.h"
87#include "MTriggerPatternDecode.h"
88#include "MBadPixelsMerge.h"
89#include "MFillH.h"
90#include "MPedCalcPedRun.h"
91#include "MPedCalcFromLoGain.h"
92#include "MFTriggerPattern.h"
93#include "MBadPixelsCalc.h"
94
95#include "MPedPhotCam.h"
96#include "MPedPhotPix.h"
97#include "MPedestalCam.h"
98#include "MPedestalPix.h"
99#include "MPedestalSubtract.h"
100
101ClassImp(MJPedestal);
102
103using namespace std;
104
105const TString MJPedestal::fgReferenceFile = "mjobs/pedestalref.rc";
106const TString MJPedestal::fgBadPixelsFile = "mjobs/badpixels_0_559.rc";
107const Float_t MJPedestal::fgExtractWinLeft = 0;
108const Float_t MJPedestal::fgExtractWinRight = 0;
109
110// --------------------------------------------------------------------------
111//
112// Default constructor.
113//
114// Sets:
115// - fExtractor to NULL,
116// - fExtractType to kUsePedRun
117// - fStorage to Normal Storage
118// - fExtractorResolution to kFALSE
119//
120MJPedestal::MJPedestal(const char *name, const char *title)
121 : fExtractor(NULL), fDisplayType(kDisplayDataCheck),
122 fExtractType(kUsePedRun), fExtractionType(kFundamental),
123 fIsUseHists(kFALSE), fDeadPixelCheck(kFALSE), fMinEvents(50)
124{
125 fName = name ? name : "MJPedestal";
126 fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
127
128 SetUsePedRun();
129 SetPathIn("");
130 SetReferenceFile();
131 SetBadPixelsFile();
132
133 SetExtractWinLeft();
134 SetExtractWinRight();
135 //
136 // Default references for case that no value references file is there
137 // (should not occur)
138 //
139
140 fPedestalMin = 4.;
141 fPedestalMax = 16.;
142 fPedRmsMin = 0.;
143 fPedRmsMax = 20.;
144 fRefPedClosedLids = 9.635;
145 fRefPedExtraGalactic = 9.93;
146 fRefPedGalactic = 10.03;
147 fRefPedRmsClosedLidsInner = 1.7;
148 fRefPedRmsExtraGalacticInner = 5.6;
149 fRefPedRmsGalacticInner = 6.92;
150 fRefPedRmsClosedLidsOuter = 1.7;
151 fRefPedRmsExtraGalacticOuter = 3.35;
152 fRefPedRmsGalacticOuter = 4.2;
153}
154
155MJPedestal::~MJPedestal()
156{
157 if (fExtractor)
158 delete fExtractor;
159}
160
161const char* MJPedestal::GetOutputFileName() const
162{
163 return Form("pedest%08d.root", fSequence.GetSequence());
164}
165
166MExtractor *MJPedestal::ReadCalibration()
167{
168 const TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence());
169
170 *fLog << inf << "Reading extractor from file: " << fname << endl;
171
172 TFile file(fname, "READ");
173 if (!file.IsOpen())
174 {
175 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
176 return NULL;
177 }
178
179 if (file.FindKey("MBadPixelsCam"))
180 {
181 MBadPixelsCam bad;
182 if (bad.Read()<=0)
183 *fLog << warn << "Unable to read MBadPixelsCam from " << fname << endl;
184 else
185 fBadPixels.Merge(bad);
186 }
187
188 if (fExtractor)
189 return fExtractor;
190
191 TObject *o=0;
192 o = file.Get("ExtractSignal");
193 if (o && !o->InheritsFrom(MExtractor::Class()))
194 {
195 *fLog << err << dbginf << "ERROR - ExtractSignal read from " << fname << " doesn't inherit from MExtractor!" << endl;
196 return NULL;
197 }
198 return o ? (MExtractor*)o->Clone("ExtractSignal") : NULL;
199}
200
201//---------------------------------------------------------------------------------
202//
203// Display the results.
204// If Display type "kDataCheck" was chosen, also the reference lines are displayed.
205//
206void MJPedestal::DisplayResult(const MParList &plist)
207{
208 if (!fDisplay)
209 return;
210
211 //
212 // Update display
213 //
214 TString title = fDisplay->GetTitle();
215 title += "-- Pedestal ";
216 title += fSequence.GetFileName();
217 title += " --";
218 fDisplay->SetTitle(title);
219
220 //
221 // Get container from list
222 //
223 MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
224 // MCalibrationPedCam &calpedcam = *(MCalibrationPedCam*)plist.FindObject("MCalibrationPedCam");
225
226 //
227 // Create container to display
228 //
229 MHCamera disp0 (geomcam, "MPedestalCam;ped", "Mean Pedestal");
230 MHCamera disp1 (geomcam, "MPedestalCam;RMS", "Pedestal RMS");
231 MHCamera disp2 (geomcam, "MCalibPedCam;histmean", "Mean Pedestal (Hist.)");
232 MHCamera disp3 (geomcam, "MCalibPedCam;histsigma", "Pedestal RMS (Hist.)");
233 MHCamera disp4 (geomcam, "MCalibPedCam;ped", "Mean Pedestal");
234 MHCamera disp5 (geomcam, "MCalibPedCam;RMS", "Pedestal RMS");
235 MHCamera disp6 (geomcam, "MCalibDiffCam;ped", "Diff. Mean Pedestal (Hist.)");
236 MHCamera disp7 (geomcam, "MCalibDiffCam;RMS", "Diff. Pedestal RMS (Hist.)");
237 MHCamera disp8 (geomcam, "MCalibDiffCam;ped", "Diff. Mean Pedestal");
238 MHCamera disp9 (geomcam, "MCalibDiffCam;AbsRMS", "Diff. Abs. Pedestal RMS");
239 MHCamera disp10(geomcam, "MCalibDiffCam;RelRMS", "Diff. Rel. Pedestal RMS");
240
241 disp0.SetCamContent(fPedestalCamOut, 0);
242 disp0.SetCamError (fPedestalCamOut, 1);
243
244 disp1.SetCamContent(fPedestalCamOut, 2);
245 disp1.SetCamError (fPedestalCamOut, 3);
246
247 /*
248 if (fIsUseHists)
249 {
250 disp2.SetCamContent(calpedcam, 0);
251 disp2.SetCamError (calpedcam, 1);
252
253 disp3.SetCamContent(calpedcam, 2);
254 disp3.SetCamError (calpedcam, 3);
255
256 disp4.SetCamContent(calpedcam, 5);
257 disp4.SetCamError (calpedcam, 6);
258
259 disp5.SetCamContent(calpedcam, 7);
260 disp5.SetCamError (calpedcam, 8);
261
262 for (UInt_t i=0;i<geomcam.GetNumPixels();i++)
263 {
264
265 MPedestalPix &ped = fPedestalCamOut[i];
266 MCalibrationPix &hist = calpedcam [i];
267 MBadPixelsPix &bad = fBadPixels[i];
268
269 if (bad.IsUnsuitable())
270 continue;
271
272 disp6.Fill(i,ped.GetPedestal()-hist.GetHiGainMean());
273 disp6.SetUsed(i);
274
275 disp7.Fill(i,hist.GetHiGainSigma()-ped.GetPedestalRms());
276 if (TMath::Abs(ped.GetPedestalRms()-hist.GetHiGainSigma()) < 4.0)
277 disp7.SetUsed(i);
278
279 disp8.Fill(i,ped.GetPedestal()-hist.GetLoGainMean());
280 disp8.SetUsed(i);
281
282 disp9.Fill(i,hist.GetLoGainSigma()-ped.GetPedestalRms());
283 if (TMath::Abs(hist.GetLoGainSigma() - ped.GetPedestalRms()) < 4.0)
284 disp9.SetUsed(i);
285 }
286 }
287 */
288
289 if (fExtractionType!=kFundamental/*fExtractorResolution*/)
290 {
291 for (UInt_t i=0;i<geomcam.GetNumPixels();i++)
292 {
293
294 MPedestalPix &pedo = fPedestalCamOut[i];
295 MPedestalPix &pedi = fPedestalCamIn[i];
296 MBadPixelsPix &bad = fBadPixels[i];
297
298 if (bad.IsUnsuitable())
299 continue;
300
301 const Float_t diff = pedo.GetPedestalRms()-pedi.GetPedestalRms();
302 const Float_t sum = 0.5*(pedo.GetPedestalRms()+pedi.GetPedestalRms());
303
304 disp9.Fill(i,pedo.GetPedestalRms()-pedi.GetPedestalRms());
305 if (pedo.IsValid() && pedi.IsValid())
306 disp9.SetUsed(i);
307
308 disp10.Fill(i,sum == 0. ? 0. : diff/sum);
309 if (pedo.IsValid() && pedi.IsValid() && sum != 0.)
310 disp10.SetUsed(i);
311 }
312 }
313
314 disp0.SetYTitle("P [cts/slice]");
315 disp1.SetYTitle("P_{rms} [cts/slice]");
316 disp2.SetYTitle("Hist. Mean [cts/slice]");
317 disp3.SetYTitle("Hist. Sigma [cts/slice]");
318 disp4.SetYTitle("Calc. Mean [cts/slice]");
319 disp5.SetYTitle("Calc. RMS [cts/slice]");
320 disp6.SetYTitle("Diff. Mean [cts/slice]");
321 disp7.SetYTitle("Diff. RMS [cts/slice]");
322 disp8.SetYTitle("Diff. Mean [cts/slice]");
323 disp9.SetYTitle("Abs.Diff.RMS [cts/slice]");
324 disp10.SetYTitle("Rel.Diff.RMS [1]");
325
326 //
327 // Display data
328 //
329 if (fDisplayType != kDisplayDataCheck && fExtractionType==kFundamental/*fExtractorResolution*/)
330 {
331 TCanvas &c3 = fDisplay->AddTab("Pedestals");
332 c3.Divide(2,3);
333
334 disp0.CamDraw(c3, 1, 2, 1);
335 disp1.CamDraw(c3, 2, 2, 6);
336 return;
337 }
338
339#if 0
340 if (fIsUseHists)
341 {
342
343 TCanvas &c3 = fDisplay->AddTab("Extractor Hist.");
344 c3.Divide(2,3);
345
346 disp2.CamDraw(c3, 1, 2, 1);
347 disp3.CamDraw(c3, 2, 2, 5);
348
349 TCanvas &c4 = fDisplay->AddTab("Extractor Calc.");
350 c4.Divide(2,3);
351
352 disp4.CamDraw(c4, 1, 2, 1);
353 disp5.CamDraw(c4, 2, 2, 5);
354
355 /*
356 TCanvas &c5 = fDisplay->AddTab("Difference Hist.");
357 c5.Divide(2,3);
358
359 disp6.CamDraw(c5, 1, 2, 1);
360 disp7.CamDraw(c5, 2, 2, 5);
361 */
362
363 TCanvas &c6 = fDisplay->AddTab("Difference Calc.");
364 c6.Divide(2,3);
365
366 disp8.CamDraw(c6, 1, 2, 1);
367 disp9.CamDraw(c6, 2, 2, 5);
368 return;
369 }
370#endif
371 if (fDisplayType == kDisplayDataCheck)
372 {
373
374 TCanvas &c3 = fDisplay->AddTab(fExtractionType!=kFundamental/*fExtractorResolution*/ ? "PedExtrd" : "Ped");
375 c3.Divide(2,3);
376
377 c3.cd(1);
378 gPad->SetBorderMode(0);
379 gPad->SetTicks();
380 MHCamera *obj1=(MHCamera*)disp0.DrawCopy("hist");
381 //
382 // for the datacheck, fix the ranges!!
383 //
384 if (fExtractionType==kFundamental/*!fExtractorResolution*/)
385 {
386 obj1->SetMinimum(fPedestalMin);
387 obj1->SetMaximum(fPedestalMax);
388 }
389 //
390 // Set the datacheck sizes:
391 //
392 FixDataCheckHist((TH1D*)obj1);
393 //
394 // set reference lines
395 //
396 DisplayReferenceLines(obj1,0);
397 //
398 // end reference lines
399 //
400 c3.cd(3);
401 gPad->SetBorderMode(0);
402 obj1->SetPrettyPalette();
403 obj1->Draw();
404
405 c3.cd(5);
406 gPad->SetBorderMode(0);
407 gPad->SetTicks();
408 TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
409 obj2->Draw();
410 obj2->SetBit(kCanDelete);
411 obj2->Fit("gaus","Q");
412 obj2->GetFunction("gaus")->SetLineColor(kYellow);
413 //
414 // Set the datacheck sizes:
415 //
416 FixDataCheckHist(obj2);
417 obj2->SetStats(1);
418
419 c3.cd(2);
420 gPad->SetBorderMode(0);
421 gPad->SetTicks();
422 MHCamera *obj3=(MHCamera*)disp1.DrawCopy("hist");
423 //
424 // for the datacheck, fix the ranges!!
425 //
426 obj3->SetMinimum(fPedRmsMin);
427 obj3->SetMaximum(fPedRmsMax);
428 //
429 // Set the datacheck sizes:
430 //
431 FixDataCheckHist((TH1D*)obj3);
432 //
433 // set reference lines
434 //
435 DisplayReferenceLines(obj3,1);
436
437 c3.cd(4);
438 gPad->SetBorderMode(0);
439 obj3->SetPrettyPalette();
440 obj3->Draw();
441
442 c3.cd(6);
443 gPad->SetBorderMode(0);
444
445 if (geomcam.InheritsFrom("MGeomCamMagic"))
446 {
447 TArrayI inner(1);
448 inner[0] = 0;
449
450 TArrayI outer(1);
451 outer[0] = 1;
452
453 TArrayI s0(6);
454 s0[0] = 6;
455 s0[1] = 1;
456 s0[2] = 2;
457 s0[3] = 3;
458 s0[4] = 4;
459 s0[5] = 5;
460
461 TArrayI s1(3);
462 s1[0] = 6;
463 s1[1] = 1;
464 s1[2] = 2;
465
466 TArrayI s2(3);
467 s2[0] = 3;
468 s2[1] = 4;
469 s2[2] = 5;
470
471 TVirtualPad *pad = gPad;
472 pad->Divide(2,1);
473
474 TH1D *inout[2];
475 inout[0] = disp1.ProjectionS(s0, inner, "Inner");
476 inout[1] = disp1.ProjectionS(s0, outer, "Outer");
477 FixDataCheckHist(inout[0]);
478 FixDataCheckHist(inout[1]);
479
480 inout[0]->SetTitle(Form("%s %s",disp1.GetTitle(),"Inner"));
481 inout[1]->SetTitle(Form("%s %s",disp1.GetTitle(),"Outer"));
482
483
484 for (int i=0; i<2; i++)
485 {
486 pad->cd(i+1);
487 gPad->SetBorderMode(0);
488 gPad->SetTicks();
489
490 inout[i]->SetDirectory(NULL);
491 inout[i]->SetLineColor(kRed+i);
492 inout[i]->SetBit(kCanDelete);
493 inout[i]->Draw();
494 inout[i]->Fit("gaus", "Q");
495
496 TLegend *leg2 = new TLegend(0.6,0.2,0.9,0.55);
497 leg2->SetHeader(inout[i]->GetName());
498 leg2->AddEntry(inout[i], inout[i]->GetName(), "l");
499
500 //
501 // Display the outliers as dead and noisy pixels
502 //
503 DisplayOutliers(inout[i]);
504
505 //
506 // Display the two half of the camera separately
507 //
508 TH1D *half[2];
509 half[0] = disp1.ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2");
510 half[1] = disp1.ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5");
511
512 for (int j=0; j<2; j++)
513 {
514 half[j]->SetLineColor(kRed+i+2*j+1);
515 half[j]->SetDirectory(NULL);
516 half[j]->SetBit(kCanDelete);
517 half[j]->Draw("same");
518 leg2->AddEntry(half[j], half[j]->GetName(), "l");
519 }
520 leg2->Draw();
521 delete leg2;
522 }
523 return;
524 }
525 }
526
527 if (fExtractionType!=kFundamental/*fExtractorResolution*/)
528 {
529
530 TCanvas &c3 = fDisplay->AddTab(fExtractionType==kWithExtractor?"PedExtrd":"PedRndm");
531 c3.Divide(2,3);
532
533 disp0.CamDraw(c3, 1, 2, 1);
534 disp1.CamDraw(c3, 2, 2, 6);
535
536 TCanvas &c13 = fDisplay->AddTab(fExtractionType==kWithExtractor?"DiffExtrd":"DiffRndm");
537 c13.Divide(2,3);
538
539 disp9.CamDraw(c13, 1, 2, 5);
540 disp10.CamDraw(c13, 2, 2, 5);
541 return;
542 }
543}
544
545
546void MJPedestal::DisplayReferenceLines(MHCamera *cam, const Int_t what) const
547{
548
549 Double_t x = cam->GetNbinsX();
550
551 const MGeomCam *geom = cam->GetGeometry();
552
553 if (geom->InheritsFrom("MGeomCamMagic"))
554 x = what ? 397 : cam->GetNbinsX();
555
556 TLine line;
557 line.SetLineStyle(kDashed);
558 line.SetLineWidth(3);
559 line.SetLineColor(kBlue);
560
561 TLegend *leg = new TLegend(0.6,0.75,0.9,0.99);
562 leg->SetBit(kCanDelete);
563
564 if (fExtractionType==kWithExtractorRndm && !(what))
565 {
566 TLine *l0 = line.DrawLine(0,0.,cam->GetNbinsX(),0.);
567 l0->SetBit(kCanDelete);
568 leg->AddEntry(l0, "Reference","l");
569 leg->Draw();
570 return;
571 }
572
573 line.SetLineColor(kBlue);
574 TLine *l1 = line.DrawLine(0, what ? fRefPedRmsGalacticInner : fRefPedGalactic,
575 x, what ? fRefPedRmsGalacticInner : fRefPedGalactic);
576 l1->SetBit(kCanDelete);
577 line.SetLineColor(kYellow);
578 TLine *l2 = line.DrawLine(0, what ? fRefPedRmsExtraGalacticInner : fRefPedExtraGalactic,
579 x, what ? fRefPedRmsExtraGalacticInner : fRefPedExtraGalactic);
580 l2->SetBit(kCanDelete);
581 line.SetLineColor(kMagenta);
582 TLine *l3 = line.DrawLine(0, what ? fRefPedRmsClosedLidsInner : fRefPedClosedLids,
583 x, what ? fRefPedRmsClosedLidsInner : fRefPedClosedLids);
584 l3->SetBit(kCanDelete);
585
586 if (geom->InheritsFrom("MGeomCamMagic"))
587 if (what)
588 {
589 const Double_t x2 = cam->GetNbinsX();
590
591 line.SetLineColor(kBlue);
592 line.DrawLine(398, fRefPedRmsGalacticOuter,
593 x2, fRefPedRmsGalacticOuter);
594
595 line.SetLineColor(kYellow);
596 line.DrawLine(398, fRefPedRmsExtraGalacticOuter,
597 x2, fRefPedRmsExtraGalacticOuter);
598
599 line.SetLineColor(kMagenta);
600 line.DrawLine(398, fRefPedRmsClosedLidsOuter,
601 x2, fRefPedRmsClosedLidsOuter);
602 }
603
604
605 leg->AddEntry(l1, "Galactic Source","l");
606 leg->AddEntry(l2, "Extra-Galactic Source","l");
607 leg->AddEntry(l3, "Closed Lids","l");
608 leg->Draw();
609}
610
611void MJPedestal::DisplayOutliers(TH1D *hist) const
612{
613 const Float_t mean = hist->GetFunction("gaus")->GetParameter(1);
614 const Float_t lolim = mean - 3.5*hist->GetFunction("gaus")->GetParameter(2);
615 const Float_t uplim = mean + 3.5*hist->GetFunction("gaus")->GetParameter(2);
616 const Stat_t dead = hist->Integral(0,hist->FindBin(lolim)-1);
617 const Stat_t noisy = hist->Integral(hist->FindBin(uplim)+1,hist->GetNbinsX()+1);
618
619 TLatex deadtex;
620 deadtex.SetTextSize(0.06);
621 deadtex.DrawLatex(0.1,hist->GetBinContent(hist->GetMaximumBin())/1.1,Form("%3i dead pixels",(Int_t)dead));
622
623 TLatex noisytex;
624 noisytex.SetTextSize(0.06);
625 noisytex.DrawLatex(0.1,hist->GetBinContent(hist->GetMaximumBin())/1.2,Form("%3i noisy pixels",(Int_t)noisy));
626}
627
628void MJPedestal::FixDataCheckHist(TH1D *hist) const
629{
630 hist->SetDirectory(NULL);
631 hist->SetStats(0);
632
633 //
634 // set the labels bigger
635 //
636 TAxis *xaxe = hist->GetXaxis();
637 TAxis *yaxe = hist->GetYaxis();
638
639 xaxe->CenterTitle();
640 yaxe->CenterTitle();
641 xaxe->SetTitleSize(0.06);
642 yaxe->SetTitleSize(0.06);
643 xaxe->SetTitleOffset(0.8);
644 yaxe->SetTitleOffset(0.5);
645 xaxe->SetLabelSize(0.05);
646 yaxe->SetLabelSize(0.05);
647}
648
649/*
650Bool_t MJPedestal::WriteEventloop(MEvtLoop &evtloop) const
651{
652 if (fOutputPath.IsNull())
653 return kTRUE;
654
655 const TString oname(GetOutputFile());
656
657 *fLog << inf << "Writing to file: " << oname << endl;
658
659 TFile file(oname, fOverwrite?"RECREATE":"NEW", "File created by MJPedestal", 9);
660 if (!file.IsOpen())
661 {
662 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
663 return kFALSE;
664 }
665
666 if (evtloop.Write(fName)<=0)
667 {
668 *fLog << err << "Unable to write MEvtloop to " << oname << endl;
669 return kFALSE;
670 }
671
672 return kTRUE;
673}
674*/
675
676void MJPedestal::SetExtractor(MExtractor* ext)
677{
678 if (ext)
679 {
680 if (fExtractor)
681 delete fExtractor;
682 fExtractor = ext ? (MExtractor*)ext->Clone(ext->GetName()) : NULL;
683 }
684 else
685 fExtractor = 0;
686}
687
688// --------------------------------------------------------------------------
689//
690// Read the following values from resource file:
691//
692// PedestalMin
693// PedestalMax
694//
695// PedRmsMin
696// PedRmsMax
697//
698// RefPedClosedLids
699// RefPedExtraGalactic
700// RefPedGalactic
701//
702// RefPedRmsClosedLidsInner
703// RefPedRmsExtraGalacticInner
704// RefPedRmsGalacticInner
705// RefPedRmsClosedLidsOuter
706// RefPedRmsExtraGalacticOuter
707// RefPedRmsGalacticOuter
708//
709void MJPedestal::ReadReferenceFile()
710{
711 TEnv refenv(fReferenceFile);
712
713 fPedestalMin = refenv.GetValue("PedestalMin",fPedestalMin);
714 fPedestalMax = refenv.GetValue("PedestalMax",fPedestalMax);
715 fPedRmsMin = refenv.GetValue("PedRmsMin",fPedRmsMin);
716 fPedRmsMax = refenv.GetValue("PedRmsMax",fPedRmsMax);
717 fRefPedClosedLids = refenv.GetValue("RefPedClosedLids",fRefPedClosedLids);
718 fRefPedExtraGalactic = refenv.GetValue("RefPedExtraGalactic",fRefPedExtraGalactic);
719 fRefPedGalactic = refenv.GetValue("RefPedGalactic",fRefPedGalactic);
720 fRefPedRmsClosedLidsInner = refenv.GetValue("RefPedRmsClosedLidsInner",fRefPedRmsClosedLidsInner);
721 fRefPedRmsExtraGalacticInner = refenv.GetValue("RefPedRmsExtraGalacticInner",fRefPedRmsExtraGalacticInner);
722 fRefPedRmsGalacticInner = refenv.GetValue("RefPedRmsGalacticInner",fRefPedRmsGalacticInner);
723 fRefPedRmsClosedLidsOuter = refenv.GetValue("RefPedRmsClosedLidsOuter",fRefPedRmsClosedLidsOuter);
724 fRefPedRmsExtraGalacticOuter = refenv.GetValue("RefPedRmsExtraGalacticOuter",fRefPedRmsExtraGalacticOuter);
725 fRefPedRmsGalacticOuter = refenv.GetValue("RefPedRmsGalacticOuter",fRefPedRmsGalacticOuter);
726}
727
728// --------------------------------------------------------------------------
729//
730// The following resource options are available:
731//
732// Do a datacheck run (read raw-data and enable display)
733// Prefix.DataCheck: Yes, No <default>
734//
735// Setup display type
736// Prefix.Display: normal <default>, datacheck, none
737//
738// Use cosmic data instead of pedestal data (DatRuns)
739// Prefix.UseData: Yes, No <default>
740//
741// Write an output file with pedestals and status-display
742// Prefix.DisableOutput: Yes, No <default>
743//
744// Name of a file containing reference values (see ReadReferenceFile)
745// Prefix.ReferenceFile: filename
746// (see ReadReferenceFile)
747//
748Bool_t MJPedestal::CheckEnvLocal()
749{
750 if (HasEnv("Display"))
751 {
752 TString type = GetEnv("Display", "normal");
753 type.ToLower();
754 if (type==(TString)"normal")
755 fDisplayType = kDisplayNormal;
756 if (type==(TString)"datacheck")
757 fDisplayType = kDisplayDataCheck;
758 if (type==(TString)"none")
759 fDisplayType = kDisplayNone;
760 }
761
762
763 SetExtractWinLeft (GetEnv("ExtractWinLeft", fExtractWinLeft ));
764 SetExtractWinRight(GetEnv("ExtractWinRight", fExtractWinRight));
765
766 fMinEvents = (UInt_t)GetEnv("MinEvents", (Int_t)fMinEvents);
767
768 if (!MJCalib::CheckEnvLocal())
769 return kFALSE;
770
771 if (HasEnv("UseData"))
772 fExtractType = GetEnv("UseData",kFALSE) ? kUseData : kUsePedRun;
773
774 if (fSequence.IsMonteCarlo() && fExtractType==kUseData)
775 {
776 // The reason is, that the standard data files contains empty
777 // (untriggered) events. If we would loop over the default 500
778 // first events of the data file you would calculate the
779 // pedestal from only some single events...
780 *fLog << inf;
781 *fLog << "Sorry, you cannot extract the starting pedestal from the first" << endl;
782 *fLog << "events in your data files... using pedestal file instead. The" << endl;
783 *fLog << "result should not differ..." << endl;
784 fExtractType = kUsePedRun;
785 }
786
787
788 if (HasEnv("UseHists"))
789 if (GetEnv("UseHists",kFALSE))
790 fIsUseHists = kTRUE;
791
792 SetNoStorage(GetEnv("DisableOutput", IsNoStorage()));
793
794 fDeadPixelCheck = GetEnv("DeadPixelsCheck", fDeadPixelCheck);
795
796 fBadPixelsFile = GetEnv("BadPixelsFile",fBadPixelsFile.Data());
797 fReferenceFile = GetEnv("ReferenceFile",fReferenceFile.Data());
798 ReadReferenceFile();
799
800 // ------------- Do not put simple resource below --------------
801
802 // Setup an environment task
803 MTaskEnv tenv("ExtractSignal");
804 tenv.SetDefault(fExtractor);
805
806 // check the resource file for it
807 if (!CheckEnv(tenv))
808 return kFALSE;
809
810// if (tenv.ReadEnv(*GetEnv(), GetEnvPrefix()+".ExtractSignal", GetEnvDebug()>2)==kERROR)
811// return kFALSE;
812
813 // If the resource file didn't change the default we are done
814 if (fExtractor==tenv.GetTask())
815 return kTRUE;
816
817 // If it changed the default check its inheritance...
818 if (!tenv.GetTask()->InheritsFrom(MExtractor::Class()))
819 {
820 *fLog << err << "ERROR: ExtractSignal from resource file doesn't inherit from MExtractor.... abort." << endl;
821 return kFALSE;
822 }
823
824 // ..and store it
825 SetExtractor((MExtractor*)tenv.GetTask());
826
827 return kTRUE;
828}
829
830//---------------------------------------------------------------------------------
831//
832Bool_t MJPedestal::WritePulsePos(TObject *obj) const
833{
834 if (IsNoStorage())
835 return kTRUE;
836
837 const TString name(Form("signal%08d.root", fSequence.GetSequence()));
838
839 TObjArray arr;
840 arr.Add(obj);
841 return WriteContainer(arr, name, fOverwrite?"RECREATE":"NEW");
842}
843
844Int_t MJPedestal::PulsePosCheck(const MParList &plist) const
845{
846 if (fIsPixelCheck)
847 {
848 MHPedestalCam *hcam = (MHPedestalCam*)plist.FindObject("MHPedestalCam");
849 if (hcam)
850 {
851 MHPedestalPix &pix1 = (MHPedestalPix&)(*hcam)[fCheckedPixId];
852 pix1.DrawClone("");
853 }
854 }
855
856 if (!fIsPulsePosCheck)
857 return kTRUE;
858
859 // FIXME:
860 // The MC cannot run over the first 2000 pedestal events since almost all
861 // events are empty, therefore a pulse pos. check is not possible, either.
862 // For the moment, have to fix the problem hardcoded...
863 //
864 // MMcEvt *evt = (MMcEvt*)plist.FindObject("MMcEvt");
865 // const Float_t meanpulsetime = evt->GetFadcTimeJitter();
866 Float_t meanpulsetime = 4.5;
867 Float_t rmspulsetime = 1.0;
868
869 MCalibrationPulseTimeCam *cam = NULL;
870 if (!fSequence.IsMonteCarlo())
871 {
872 if (fIsPixelCheck)
873 {
874 MHCalibrationPulseTimeCam *hcam = (MHCalibrationPulseTimeCam*)plist.FindObject("MHCalibrationPulseTimeCam");
875 if (!hcam)
876 {
877 *fLog << err << "MHCalibrationPulseTimeCam not found... abort." << endl;
878 return kFALSE;
879 }
880 hcam->DrawClone();
881 gPad->SaveAs(Form("%s/PulsePosTest_all.root",fPathOut.Data()));
882
883 MHCalibrationPix &pix = (*hcam)[fCheckedPixId];
884 pix.DrawClone();
885 gPad->SaveAs(Form("%s/PulsePosTest_Pixel%04d.root",fPathOut.Data(),fCheckedPixId));
886 }
887
888 cam = (MCalibrationPulseTimeCam*)plist.FindObject("MCalibrationPulseTimeCam");
889 if (!cam)
890 {
891 *fLog << err << "MCalibrationPulseTimeCam not found... abort." << endl;
892 return kFALSE;
893 }
894
895 meanpulsetime = cam->GetAverageArea(0).GetHiGainMean();
896 rmspulsetime = cam->GetAverageArea(0).GetHiGainRms();
897 }
898
899 if (!WritePulsePos(cam))
900 return kFALSE;
901
902 *fLog << all << "Mean pulse time/Avg pos.of maximum (" << (fSequence.IsMonteCarlo()?"MC":"cosmics") << "): ";
903 *fLog << meanpulsetime << "+-" << rmspulsetime << endl;
904
905 MExtractTimeAndCharge *ext = dynamic_cast<MExtractTimeAndCharge*>(fExtractor);
906 if (!ext)
907 {
908 *fLog << warn << "WARNING - no extractor found inheriting from MExtractTimeAndCharge... no pulse position check." << endl;
909 return kTRUE;
910 }
911
912 const Int_t hi0 = ext->GetHiGainFirst();
913 const Int_t lo1 = ext->GetLoGainLast();
914 Int_t hi1 = ext->GetHiGainLast();
915 Int_t lo0 = ext->GetLoGainFirst();
916
917 //
918 // This is for data without lo-gains
919 //
920 const Bool_t haslo = ext->HasLoGain();
921
922 //
923 // Get the ranges for the new extractor setting. The window
924 // size is always rounded to the next higher integer.
925 //
926 const Int_t wshigain = ext->GetWindowSizeHiGain();
927 const Int_t wslogain = ext->GetWindowSizeLoGain();
928
929 //
930 // Here we calculate the end of the lo-gain range
931 // as it is done in MExtractTimeAndCharge
932 //
933 const Double_t poshi = meanpulsetime;
934 const Double_t poslo = poshi + ext->GetOffsetLoGain();
935 const Double_t poslo2 = poslo + ext->GetLoGainStartShift();
936
937 //
938 // Do the right side checks range checks
939 //
940 if (poshi+wshigain+fExtractWinRight > hi1-0.5)
941 {
942 *fLog << err;
943 *fLog << "ERROR - Pulse is too much to the right, out of hi-gain range [";
944 *fLog << hi0 << "," << hi1 << "]" << endl;
945 *fLog << endl;
946 return -2;
947 }
948
949 if (haslo && poslo+wslogain+fExtractWinRight > lo1-0.5)
950 {
951 *fLog << err;
952 *fLog << "ERROR - Pulse is too much to the right, out of lo-gain range [";
953 *fLog << lo0 << "," << lo1 << "]" << endl;
954 return -2;
955 }
956
957 //
958 // Do the left side checks range checks
959 //
960 if (poshi-fExtractWinLeft < hi0+0.5)
961 {
962 *fLog << err;
963 *fLog << "ERROR - Pulse is too much to the left, out of hi-gain range [";
964 *fLog << hi0 << "," << hi1 << "]" << endl;
965 return -3;
966 }
967
968 if (haslo && poslo2-fExtractWinLeft < lo0+0.5)
969 {
970 *fLog << warn;
971 *fLog << "WARNING - Pulse is too much to the left, out of lo-gain range [";
972 *fLog << lo0 << "," << lo1 << "]" << endl;
973 *fLog << "Trying to match extraction window and pulse position..." << endl;
974
975 //
976 // Set and store the new ranges
977 //
978 Int_t shift = 0;
979 while (poslo2-fExtractWinLeft < lo0+0.5)
980 {
981 hi1--;
982 lo0--;
983
984 if (poshi+wshigain+fExtractWinRight > hi1-0.5)
985 {
986 *fLog << err << "ERROR - No proper extraction window found.... abort." << endl;
987 return -3;
988 }
989 }
990
991 *fLog << "Changed extraction window by " << shift;
992 *fLog << "slices to hi-gain [" << hi0 << "," << hi1;
993 *fLog << "] and lo-gain [" << lo0 << "," << lo1 << "]" << endl;
994
995 ext->SetRange(hi0, hi1, lo0, lo1);
996 }
997
998 return kTRUE;
999}
1000
1001Int_t MJPedestal::Process()
1002{
1003 if (!fSequence.IsValid())
1004 {
1005 *fLog << err << "ERROR - Sequence invalid..." << endl;
1006 return kFALSE;
1007 }
1008
1009 // --------------------------------------------------------------------------------
1010
1011 const TString type = IsUseData() ? "data" : "pedestal";
1012
1013 *fLog << inf;
1014 fLog->Separator(GetDescriptor());
1015 *fLog << "Calculate MPedestalCam from " << type << "-runs ";
1016 *fLog << fSequence.GetFileName() << endl;
1017 *fLog << endl;
1018
1019 // --------------------------------------------------------------------------------
1020
1021 if (!CheckEnv())
1022 return kFALSE;
1023
1024 MParList plist;
1025 MTaskList tlist;
1026 plist.AddToList(&tlist);
1027 plist.AddToList(this); // take care of fDisplay!
1028
1029 MReadMarsFile read("Events");
1030 MRawFileRead rawread(NULL);
1031
1032 MDirIter iter;
1033 if (fSequence.IsValid())
1034 {
1035 const Int_t n0 = IsUseData()
1036 ? fSequence.SetupDatRuns(iter, 0, !fSequence.IsMonteCarlo())
1037 : fSequence.SetupPedRuns(iter, 0, !fSequence.IsMonteCarlo());
1038
1039 if (n0<=0)
1040 return kFALSE;
1041 }
1042
1043 if (!fSequence.IsMonteCarlo())
1044 {
1045 rawread.AddFiles(iter);
1046 tlist.AddToList(&rawread);
1047 }
1048 else
1049 {
1050 read.DisableAutoScheme();
1051 read.AddFiles(iter);
1052 tlist.AddToList(&read);
1053 }
1054
1055 // Setup Tasklist
1056 plist.AddToList(&fPedestalCamOut);
1057 plist.AddToList(&fBadPixels);
1058
1059 MGeomApply geomapl;
1060 MBadPixelsMerge merge(&fBadPixels);
1061
1062 MPedCalcPedRun pedcalc;
1063 //pedcalc.SetPedestalUpdate(kFALSE);
1064
1065 MPedCalcFromLoGain pedlogain;
1066 pedlogain.SetPedestalUpdate(kFALSE);
1067
1068 MHPedestalCam hpedcam;
1069 hpedcam.SetPedestalsOut(&fPedestalCamOut);
1070 if (fExtractionType != kFundamental)
1071 hpedcam.SetRenorm(kTRUE);
1072
1073 // To have it in the parlist for MEnv!
1074 MHCalibrationPulseTimeCam pulcam;
1075 plist.AddToList(&pulcam);
1076 MFillH fillpul(&pulcam, "MPedestalSubtractedEvt", "FillPulseTime");
1077 fillpul.SetBit(MFillH::kDoNotDisplay);
1078
1079 tlist.AddToList(&geomapl);
1080 tlist.AddToList(&merge);
1081
1082 if (!fPathIn.IsNull())
1083 {
1084 fExtractor = ReadCalibration();
1085 if (!fExtractor)
1086 return kFALSE;
1087
1088 // The requested setup might have been overwritten
1089 if (!CheckEnv(*fExtractor))
1090 return kFALSE;
1091
1092 *fLog << all;
1093 *fLog << underline << "Signal Extractor found in calibration file and setup:" << endl;
1094 fExtractor->Print();
1095 *fLog << endl;
1096 }
1097
1098 //
1099 // Read bad pixels from outside
1100 //
1101 if (!fBadPixelsFile.IsNull())
1102 {
1103 *fLog << inf << "Excluding: " << fBadPixelsFile << endl;
1104 ifstream fin(fBadPixelsFile);
1105 fBadPixels.AsciiRead(fin);
1106 }
1107
1108 // This will make that for data with version less than 5, where
1109 // trigger patterns were not yet correct, all the events in the real
1110 // data file will be processed. In any case there are no interleaved
1111 // calibration events in such data, so this is fine.
1112 MTriggerPatternDecode decode;
1113 MFTriggerPattern fcalib("CalibFilter");
1114 fcalib.SetDefault(kTRUE);
1115 fcalib.DenyCalibration();
1116 fcalib.DenyPedestal();
1117
1118 tlist.AddToList(&decode);
1119
1120 // produce pedestal subtracted raw-data
1121 MPedestalSubtract pedsub;
1122 if (fExtractor && fExtractionType!=kFundamental)
1123 pedsub.SetPedestalCam(&fPedestalCamIn);
1124 else
1125 pedsub.SetNamePedestalCam(""); // Only copy hi- and lo-gain together!
1126 tlist.AddToList(&pedsub);
1127
1128 if (fIsPulsePosCheck)
1129 {
1130 fillpul.SetFilter(&fcalib);
1131 tlist.AddToList(&fcalib);
1132 tlist.AddToList(&fillpul);
1133 }
1134
1135 //MFillH fillC("MHPedestalCor", "MPedestalSubtractedEvt", "FillAcor");
1136 //fillC.SetNameTab("Acor");
1137 //if (fExtractor && fExtractionType==kWithExtractorRndm)
1138 // tlist.AddToList(&fillC);
1139
1140 // ----------------------------------------------------------------------
1141 // Now we make sure, that in all cases the ranges are setup correctly
1142 // ----------------------------------------------------------------------
1143 MTaskEnv taskenv("ExtractPedestal");
1144
1145 //------------------------------
1146 MFTriggerPattern ftp2("PedestalFilter");
1147 ftp2.SetDefault(kTRUE);
1148 ftp2.DenyCalibration();
1149 ftp2.RequirePedestal();
1150
1151 if (!fSequence.IsMonteCarlo() && (!fExtractor || !fExtractor->HasLoGain()))
1152 {
1153 taskenv.SetFilter(&ftp2);
1154 tlist.AddToList(&ftp2);
1155 }
1156 //------------------------------
1157
1158 switch (fExtractType)
1159 {
1160 case kUsePedRun:
1161 // In case other than 'fundamental' second argument is obsolete
1162 // pedcalc.SetExtractWindow(0,14); // kUsePedRun (take default from class)
1163 taskenv.SetDefault(&pedcalc);
1164 tlist.AddToList(&taskenv);
1165 break;
1166
1167 case kUseData:
1168 // In case other than 'fundamental' second argument is obsolete
1169 // pedlogain.SetExtractWindow(15,14); // kUseData (take default from class)
1170 taskenv.SetDefault(&pedlogain);
1171 tlist.AddToList(&taskenv);
1172 break;
1173 }
1174
1175 if (fIsUseHists && fExtractor)
1176 {
1177 if (fExtractor->InheritsFrom("MExtractTimeAndCharge"))
1178 {
1179 if (fExtractionType!=kFundamental)
1180 {
1181 const MExtractTimeAndCharge &e = *static_cast<MExtractTimeAndCharge*>(fExtractor);
1182 hpedcam.SetFitStart(-5*e.GetWindowSizeHiGain());
1183 }
1184 else
1185 hpedcam.SetFitStart(10.);
1186 }
1187
1188 plist.AddToList(&hpedcam);
1189 }
1190
1191 pedcalc.SetPedestalsOut(&fPedestalCamOut);
1192 pedlogain.SetPedestalsOut(&fPedestalCamOut);
1193
1194 // kFundamental
1195 if (fExtractor)
1196 {
1197 if (fExtractionType!=kFundamental)
1198 {
1199 pedcalc.SetRandomCalculation(fExtractionType==kWithExtractorRndm);
1200 pedlogain.SetRandomCalculation(fExtractionType==kWithExtractorRndm);
1201
1202 pedcalc.SetExtractor((MExtractTimeAndCharge*)fExtractor);
1203 pedlogain.SetExtractor((MExtractTimeAndCharge*)fExtractor);
1204 }
1205 else
1206 {
1207 pedcalc.SetRangeFromExtractor(*fExtractor);
1208 pedlogain.SetRangeFromExtractor(*fExtractor);
1209 }
1210
1211 if (!fExtractor->InheritsFrom("MExtractTimeAndCharge") && fExtractionType!=kFundamental)
1212 {
1213 *fLog << inf;
1214 *fLog << "Signal extractor doesn't inherit from MExtractTimeAndCharge..." << endl;
1215 *fLog << " --> falling back to fundamental pedestal extraction." << endl;
1216 fExtractionType=kFundamental;
1217 }
1218 }
1219 else
1220 {
1221 *fLog << warn << GetDescriptor() << ": WARNING - No extractor has been handed over! " << endl;
1222 *fLog << "Taking default window for pedestal extraction. The calculated pedestal RMS" << endl;
1223 *fLog << "will probably not match with future pedestal RMS' from different extraction" << endl;
1224 *fLog << "windows." << endl;
1225 }
1226
1227
1228 /*
1229 MHCamEvent evt0(0, "Ped", "Pedestal;;P [cnts/sl]");
1230 MHCamEvent evt1(2, "PedRms", "Pedestal RMS;;\\sigma_{p} [cnts/sl]");
1231
1232 MFillH fill0(&evt0, &fPedestalCamOut, "FillPedestal");
1233 MFillH fill1(&evt1, &fPedestalCamOut, "FillPedRms");
1234
1235 tlist.AddToList(&fill0);
1236 tlist.AddToList(&fill1);
1237 */
1238
1239 //
1240 // Execute the eventloop
1241 //
1242 MEvtLoop evtloop(fName);
1243 evtloop.SetParList(&plist);
1244 evtloop.SetDisplay(fDisplay);
1245 evtloop.SetLogStream(fLog);
1246 if (!SetupEnv(evtloop))
1247 return kFALSE;
1248
1249 // if (!WriteEventloop(evtloop))
1250 // return kFALSE;
1251
1252 // Execute first analysis
1253 if (!evtloop.Eventloop(fMaxEvents))
1254 {
1255 *fLog << err << GetDescriptor() << ": Failed." << endl;
1256 return kFALSE;
1257 }
1258
1259 if (taskenv.GetNumExecutions()<fMinEvents)
1260 {
1261 *fLog << err << GetDescriptor() << ": Failed. Less than the required " << fMinEvents << " evts processed." << endl;
1262 return kFALSE;
1263 }
1264
1265 if (fDeadPixelCheck)
1266 {
1267 MBadPixelsCalc calc;
1268 calc.SetPedestalLevelVarianceLo(4.5);
1269 calc.SetPedestalLevelVarianceHi();
1270 calc.SetPedestalLevel();
1271 if (!CheckEnv(calc))
1272 return kFALSE;
1273 calc.SetGeomCam(dynamic_cast<MGeomCam*>(plist.FindObject("MGeomCam")));
1274 if (!calc.CheckPedestalRms(fBadPixels, fPedestalCamOut))
1275 {
1276 *fLog << err << "ERROR - MBadPixelsCalc::CheckPedestalRms failed...." << endl;
1277 return kFALSE;
1278 }
1279 }
1280
1281 if (fDisplayType!=kDisplayNone)
1282 DisplayResult(plist);
1283
1284 // if (!WriteResult())
1285 // return kFALSE;
1286
1287 const Int_t rc = PulsePosCheck(plist);
1288 if (rc<1)
1289 return rc;
1290
1291 *fLog << all << GetDescriptor() << ": Done." << endl;
1292 *fLog << endl << endl;
1293
1294 return kTRUE;
1295}
Note: See TracBrowser for help on using the repository browser.