source: trunk/Mars/mhist/MHCamera.cc@ 19345

Last change on this file since 19345 was 19345, checked in by tbretz, 3 years ago
Improves the behaviour with RecursiveRemove. Strictly speaking this change might only be necessary if a class contains more than one member which is bound to recursive remove. Onth other hand setting all members to NULL which might be affected by RecursiveRemove is not wrong either.
File size: 73.5 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, 05/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Harald Kornmayer, 1/2001
20! Author(s): Markus Gaug, 03/2004 <mailto:markus@ifae.es>
21!
22! Copyright: MAGIC Software Development, 2000-2008
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// MHCamera
30//
31// Camera Display, based on a TH1D. Pleas be carefull using the
32// underlaying TH1D.
33//
34// To change the scale to a logarithmic scale SetLogy() of the Pad.
35//
36// You can correct for the abberation. Assuming that the distance
37// between the mean position of the light distribution and the position
38// of a perfect reflection on a perfect mirror in the distance r on
39// the camera plane is dr it is d = a*dr while a is the abberation
40// constant (for the MAGIC mirror it is roughly 0.0713). You can
41// set this constant by calling SetAbberation(a) which results in a
42// 'corrected' display (all outer pixels are shifted towards the center
43// of the camera to correct for this abberation)
44//
45// Be carefull: Entries in this context means Entries/bin or Events
46//
47// FIXME? Maybe MHCamera can take the fLog object from MGeomCam?
48//
49////////////////////////////////////////////////////////////////////////////
50#include "MHCamera.h"
51
52#include <fstream>
53#include <iostream>
54
55#include <TBox.h>
56#include <TArrow.h>
57#include <TLatex.h>
58#include <TStyle.h>
59#include <TCanvas.h>
60#include <TArrayF.h>
61#include <TRandom.h>
62#include <TPaveText.h>
63#include <TPaveStats.h>
64#include <TClonesArray.h>
65#include <THistPainter.h>
66#include <THLimitsFinder.h>
67#include <TProfile.h>
68#include <TH1.h>
69#include <TF1.h>
70#include <TCanvas.h>
71#include <TLegend.h>
72
73#include "MLog.h"
74#include "MLogManip.h"
75
76#include "MH.h"
77#include "MString.h"
78#include "MBinning.h"
79
80#include "MGeom.h"
81#include "MGeomCam.h"
82
83#include "MCamEvent.h"
84
85#include "MArrayD.h"
86#include "MMath.h" // MMath::GaussProb
87
88ClassImp(MHCamera);
89
90using namespace std;
91
92void MHCamera::Init()
93{
94 Sumw2();
95
96 UseCurrentStyle();
97
98 SetDirectory(NULL);
99
100 SetLineColor(kGreen);
101 SetMarkerStyle(kFullDotMedium);
102 SetXTitle("Pixel Index");
103
104 fNotify = new TList;
105 fNotify->SetBit(kMustCleanup);
106 gROOT->GetListOfCleanups()->Add(fNotify);
107
108 /*
109 TVirtualPad *save = gPad;
110 gPad = 0;
111#if ROOT_VERSION_CODE < ROOT_VERSION(3,01,06)
112 SetPalette(1, 0);
113#endif
114 */
115/*
116#if ROOT_VERSION_CODE < ROOT_VERSION(4,04,00)
117 SetPrettyPalette();
118#elese
119 // WORAROUND - FIXME: Calling it many times becomes slower and slower
120 SetInvDeepBlueSeaPalette();
121#endif
122 gPad = save;
123*/
124}
125
126// ------------------------------------------------------------------------
127//
128// Default Constructor. To be used by the root system ONLY.
129//
130MHCamera::MHCamera() : TH1D(), fGeomCam(NULL), fAbberation(0)
131{
132 Init();
133}
134
135// ------------------------------------------------------------------------
136//
137// Constructor. Makes a clone of MGeomCam. Removed the TH1D from the
138// current directory. Calls Sumw2(). Set the histogram line color
139// (for error bars) to Green and the marker style to kFullDotMedium.
140//
141MHCamera::MHCamera(const MGeomCam &geom, const char *name, const char *title)
142: fGeomCam(NULL), fAbberation(0)
143{
144 //fGeomCam = (MGeomCam*)geom.Clone();
145 SetGeometry(geom, name, title);
146 Init();
147
148 //
149 // root 3.02
150 // * base/inc/TObject.h:
151 // register BIT(8) as kNoContextMenu. If an object has this bit set it will
152 // not get an automatic context menu when clicked with the right mouse button.
153}
154
155// ------------------------------------------------------------------------
156//
157// Clone the MHCamera via TH1D::Clone and make sure that the new object is
158// not removed from the current directory.
159//
160TObject *MHCamera::Clone(const char *newname) const
161{
162 MHCamera *rc = static_cast<MHCamera*>(TH1D::Clone(newname));
163 rc->SetDirectory(NULL);
164
165 // fGeomCam need special treatment due to its TObjArray
166 if (rc->fGeomCam && fGeomCam)
167 {
168 delete rc->fGeomCam;
169 rc->fGeomCam = static_cast<MGeomCam*>(fGeomCam->Clone());
170 }
171
172 return rc;
173}
174
175void MHCamera::SetGeometry(const MGeomCam &geom, const char *name, const char *title)
176{
177 SetNameTitle(name, title);
178
179 TAxis &x = *GetXaxis();
180
181 SetBins(geom.GetNumPixels(), 0, 1);
182 x.Set(geom.GetNumPixels(), -0.5, geom.GetNumPixels()-0.5);
183
184 //SetBins(geom.GetNumPixels(), -0.5, geom.GetNumPixels()-0.5);
185 //Rebuild();
186
187 if (fGeomCam)
188 delete fGeomCam;
189 fGeomCam = (MGeomCam*)geom.Clone();
190
191 fUsed.Set(geom.GetNumPixels());
192 fUsed.Reset();
193
194 fBinEntries.Set(geom.GetNumPixels()+2);
195 fBinEntries.Reset();
196}
197
198// ------------------------------------------------------------------------
199//
200// Destructor. Deletes the cloned fGeomCam and the notification list.
201//
202MHCamera::~MHCamera()
203{
204 if (fGeomCam)
205 {
206 delete fGeomCam;
207 fGeomCam = 0;
208 }
209 if (fNotify)
210 {
211 delete fNotify;
212 fNotify = 0;
213 }
214}
215
216// ------------------------------------------------------------------------
217//
218// Return kTRUE for sector<0. Otherwise return kTRUE only if the specified
219// sector idx matches the sector of the pixel with index idx.
220//
221Bool_t MHCamera::MatchSector(Int_t idx, const TArrayI &sector, const TArrayI &aidx) const
222{
223 const MGeom &pix = (*fGeomCam)[idx];
224 return FindVal(sector, pix.GetSector()) && FindVal(aidx, pix.GetAidx());
225}
226
227// ------------------------------------------------------------------------
228//
229// Taken from TH1D::Fill(). Uses the argument directly as bin index.
230// Doesn't increment the number of entries.
231//
232// -*-*-*-*-*-*-*-*Increment bin with abscissa X by 1*-*-*-*-*-*-*-*-*-*-*
233// ==================================
234//
235// if x is less than the low-edge of the first bin, the Underflow bin is incremented
236// if x is greater than the upper edge of last bin, the Overflow bin is incremented
237//
238// If the storage of the sum of squares of weights has been triggered,
239// via the function Sumw2, then the sum of the squares of weights is incremented
240// by 1 in the bin corresponding to x.
241//
242// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
243Int_t MHCamera::Fill(Axis_t x)
244{
245#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
246 if (fBuffer) return BufferFill(x,1);
247#endif
248 const Int_t bin = (Int_t)x+1;
249 AddBinContent(bin);
250 fBinEntries[bin]++;
251 if (fSumw2.fN)
252 fSumw2.fArray[bin]++;
253
254 if (bin<=0 || bin>fNcells-2)
255 return -1;
256
257 fTsumw++;
258 fTsumw2++;
259 fTsumwx += x;
260 fTsumwx2 += x*x;
261 return bin;
262}
263
264// ------------------------------------------------------------------------
265//
266// Taken from TH1D::Fill(). Uses the argument directly as bin index.
267// Doesn't increment the number of entries.
268//
269// -*-*-*-*-*-*Increment bin with abscissa X with a weight w*-*-*-*-*-*-*-*
270// =============================================
271//
272// if x is less than the low-edge of the first bin, the Underflow bin is incremented
273// if x is greater than the upper edge of last bin, the Overflow bin is incremented
274//
275// If the storage of the sum of squares of weights has been triggered,
276// via the function Sumw2, then the sum of the squares of weights is incremented
277// by w^2 in the bin corresponding to x.
278//
279// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
280Int_t MHCamera::Fill(Axis_t x, Stat_t w)
281{
282#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
283 if (fBuffer) return BufferFill(x,w);
284#endif
285 const Int_t bin = (Int_t)x+1;
286 if (bin<0 || bin>fNcells-1)
287 return -1;
288
289 AddBinContent(bin, w);
290 fBinEntries[bin]++;
291 if (fSumw2.fN)
292 fSumw2.fArray[bin] += w*w;
293
294 if (bin<=0 || bin>fNcells-2)
295 return -1;
296
297 const Stat_t z = (w > 0 ? w : -w);
298 fTsumw += z;
299 fTsumw2 += z*z;
300 fTsumwx += z*x;
301 fTsumwx2 += z*x*x;
302 return bin;
303}
304
305// ------------------------------------------------------------------------
306//
307// Use x and y in millimeters
308//
309Int_t MHCamera::Fill(Axis_t x, Axis_t y, Stat_t w)
310{
311 if (fNcells<=1 || IsFreezed())
312 return -1;
313
314 for (Int_t idx=0; idx<fNcells-2; idx++)
315 {
316 if (!(*fGeomCam)[idx].IsInside(x, y))
317 continue;
318
319 SetUsed(idx);
320 return Fill(idx, w);
321 }
322 return -1;
323}
324
325// ------------------------------------------------------------------------
326//
327// Call this if you want to change the display status (displayed or not)
328// for all pixels. val==0 means that the pixel is not displayed.
329//
330void MHCamera::SetUsed(const TArrayC &arr)
331{
332 if (fNcells-2 != arr.GetSize())
333 {
334 gLog << warn << "WARNING - MHCamera::SetUsed: array size mismatch... ignored." << endl;
335 return;
336 }
337
338 for (Int_t idx=0; idx<fNcells-2; idx++)
339 arr[idx] ? SetUsed(idx) : ResetUsed(idx);
340}
341
342// ------------------------------------------------------------------------
343//
344// Return the mean value of all entries which are used if all=kFALSE and
345// of all entries if all=kTRUE if sector<0. If sector>=0 only
346// entries with match the given sector are taken into account.
347//
348Stat_t MHCamera::GetMeanSectors(const TArrayI &sector, const TArrayI &aidx, Bool_t ball) const
349{
350 if (fNcells<=1)
351 return 0;
352
353 Int_t n=0;
354
355 Stat_t mean = 0;
356 for (int i=0; i<fNcells-2; i++)
357 {
358 if ((ball || IsUsed(i)) && MatchSector(i, sector, aidx))
359 {
360 if (TestBit(kProfile) && fBinEntries[i+1]==0)
361 continue;
362
363 mean += TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
364 n++;
365 }
366 }
367
368 return n==0 ? 0 : mean/n;
369}
370
371UInt_t MHCamera::GetNumUsedSector(const TArrayI &sector, const TArrayI &aidx) const
372{
373 if (fNcells<=1)
374 return 0;
375
376 Int_t n=0;
377
378 for (int i=0; i<fNcells-2; i++)
379 {
380 if (!IsUsed(i) || !MatchSector(i, sector, aidx))
381 continue;
382
383 if (TestBit(kProfile) && fBinEntries[i+1]==0)
384 continue;
385 n++;
386 }
387
388 // return Median of the profile data
389 return n;
390}
391
392// ------------------------------------------------------------------------
393//
394// Return the median value of all entries which are used if all=kFALSE and
395// of all entries if all=kTRUE if sector<0. If sector>=0 only
396// entries with match the given sector are taken into account.
397//
398Stat_t MHCamera::GetMedianSectors(const TArrayI &sector, const TArrayI &aidx, Bool_t ball) const
399{
400 if (fNcells<=1)
401 return 0;
402
403 TArrayD arr(fNcells-2);
404 Int_t n=0;
405
406 for (int i=0; i<fNcells-2; i++)
407 {
408 if ((ball || IsUsed(i)) && MatchSector(i, sector, aidx))
409 {
410 if (TestBit(kProfile) && fBinEntries[i+1]==0)
411 continue;
412
413 arr[n++] = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
414 }
415 }
416
417 // return Median of the profile data
418 return TMath::Median(n, arr.GetArray());
419}
420
421// ------------------------------------------------------------------------
422//
423// Return the sqrt variance of all entries which are used if all=kFALSE and
424// of all entries if all=kTRUE if sector<0. If sector>=0 only
425// entries with match the given sector are taken into account.
426//
427Stat_t MHCamera::GetRmsSectors(const TArrayI &sector, const TArrayI &aidx, Bool_t ball) const
428{
429 if (fNcells<=1)
430 return -1;
431
432 Int_t n=0;
433
434 Stat_t sum = 0;
435 Stat_t sq = 0;
436 for (int i=0; i<fNcells-2; i++)
437 {
438 if ((ball || IsUsed(i)) && MatchSector(i, sector, aidx))
439 {
440 if (TestBit(kProfile) && fBinEntries[i+1]==0)
441 continue;
442
443 const Double_t val = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
444
445 sum += val;
446 sq += val*val;
447 n++;
448 }
449 }
450
451 if (n==0)
452 return 0;
453
454 sum /= n;
455 sq /= n;
456
457 return TMath::Sqrt(sq-sum*sum);
458}
459
460// ------------------------------------------------------------------------
461//
462// Return the median value (divided by MMath::GausProb(1.0)=68.3%) of the
463// distribution of abs(y[i]-Median). This is my Median equivalent of the RMS.
464// Return the deviation of all entries which are used if all=kFALSE and
465// of all entries if all=kTRUE if sector<0. If sector>=0 only
466// entries with match the given sector are taken into account.
467//
468Stat_t MHCamera::GetDevSectors(const TArrayI &sector, const TArrayI &aidx, Bool_t ball) const
469{
470 if (fNcells<=1)
471 return 0;
472
473 TArrayD arr(fNcells-2);
474 Int_t n=0;
475
476 for (int i=0; i<fNcells-2; i++)
477 {
478 if ((ball || IsUsed(i)) && MatchSector(i, sector, aidx))
479 {
480 if (TestBit(kProfile) && fBinEntries[i+1]==0)
481 continue;
482
483 arr[n++] = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
484 }
485 }
486
487 // return Median of the profile data
488 return MMath::MedianDev(n, arr.GetArray());
489}
490
491// ------------------------------------------------------------------------
492//
493// Return the minimum contents of all pixels (if all is set, otherwise
494// only of all 'used' pixels), fMinimum if fMinimum set. If sector>=0
495// only pixels with matching sector number are taken into account.
496//
497Double_t MHCamera::GetMinimumSectors(const TArrayI &sector, const TArrayI &aidx, Bool_t ball) const
498{
499 if (fMinimum != -1111)
500 return fMinimum;
501
502 if (fNcells<=1)
503 return 0;
504
505 Double_t minimum=FLT_MAX;
506
507 for (Int_t i=0; i<fNcells-2; i++)
508 {
509 if (TestBit(kProfile) && fBinEntries[i+1]==0)
510 continue;
511
512 const Double_t val = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
513 if (MatchSector(i, sector, aidx) && (ball || IsUsed(i)) && val<minimum && TMath::Finite(val))
514 minimum = val;
515 }
516
517 return minimum;
518}
519
520// ------------------------------------------------------------------------
521//
522// Return the maximum contents of all pixels (if all is set, otherwise
523// only of all 'used' pixels), fMaximum if fMaximum set. If sector>=0
524// only pixels with matching sector number are taken into account.
525//
526Double_t MHCamera::GetMaximumSectors(const TArrayI &sector, const TArrayI &aidx, Bool_t ball) const
527{
528 if (fMaximum!=-1111)
529 return fMaximum;
530
531 if (fNcells<=1)
532 return 1;
533
534 Double_t maximum=-FLT_MAX;
535 for (Int_t i=0; i<fNcells-2; i++)
536 {
537 if (TestBit(kProfile) && fBinEntries[i+1]==0)
538 continue;
539
540 const Double_t val = TestBit(kProfile) ? fArray[i+1]/fBinEntries[i+1] : fArray[i+1];
541 if (MatchSector(i, sector, aidx) && (ball || IsUsed(i)) && val>maximum && TMath::Finite(val))
542 maximum = val;
543 }
544
545 return maximum;
546}
547
548// ------------------------------------------------------------------------
549//
550// Return the number of bins (excluding under- and overflow) for which
551// GetBinContent returns a value > t
552//
553Int_t MHCamera::GetNumBinsAboveThreshold(Double_t t) const
554{
555 Int_t n = 0;
556 for (Int_t i=0; i<fNcells-2; i++)
557 if (GetBinContent(i+1)>t)
558 n++;
559
560 return n;
561}
562
563// ------------------------------------------------------------------------
564//
565// Return the number of bins (excluding under- and overflow) for which
566// GetBinContent returns a value < t
567//
568Int_t MHCamera::GetNumBinsBelowThreshold(Double_t t) const
569{
570 Int_t n = 0;
571 for (Int_t i=0; i<fNcells-2; i++)
572 if (GetBinContent(i+1)>t)
573 n++;
574
575 return n;
576}
577
578// ------------------------------------------------------------------------
579//
580// Call this function to draw the camera layout into your canvas.
581// Setup a drawing canvas. Add this object and all child objects
582// (hexagons, etc) to the current pad. If no pad exists a new one is
583// created. (To access the 'real' pad containing the camera you have
584// to do a cd(1) in the current layer.
585//
586// To draw a camera into its own pad do something like:
587//
588// MGeomCamMagic m;
589// MHCamera *d=new MHCamera(m);
590//
591// TCanvas *c = new TCanvas;
592// c->Divide(2,1);
593// c->cd(1);
594//
595// d->FillRandom();
596// d->Draw();
597// d->SetBit(kCanDelete);
598//
599// There are several drawing options:
600// 'hist' Draw as a standard TH1 histogram (value vs. pixel index)
601// 'box' Draw hexagons which size is in respect to its contents
602// 'nocol' Leave the 'boxed' hexagons empty
603// 'pixelindex' Display the pixel index in each pixel
604// 'sectorindex' Display the sector index in each pixel
605// 'content' Display the relative content aligned to GetMaximum() and
606// GeMinimum() ((val-min)/(max-min))
607// 'text' Draw GetBinContent as char
608// 'integer' Draw GetBinContent as integer (TMath::Nint)
609// 'float' Draw GetBinContent as float
610// 'proj' Display the y-projection of the histogram
611// 'pal0' Use Pretty palette
612// 'pal1' Use Deep Blue Sea palette
613// 'pal2' Use Inverse Depp Blue Sea palette
614// 'same' Draw trandparent pixels on top of an existing pad. This
615// makes it possible to draw the camera image on top of an
616// existing TH2, but also allows for distorted camera images
617//
618void MHCamera::Draw(Option_t *option)
619{
620 const Bool_t hassame = TString(option).Contains("same", TString::kIgnoreCase) && gPad;
621
622 // root 3.02:
623 // gPad->SetFixedAspectRatio()
624 const Color_t col = gPad ? gPad->GetFillColor() : 16;
625 TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas("CamDisplay", "Mars Camera Display", 656, 600);
626
627 if (!hassame)
628 {
629 pad->SetBorderMode(0);
630 pad->SetFillColor(col);
631
632 //
633 // Create an own pad for the MHCamera-Object which can be
634 // resized in paint to keep the correct aspect ratio
635 //
636 // The margin != 0 is a workaround for a problem in root 4.02/00
637 //pad->Divide(1, 1, 1e-10, 1e-10, col);
638 //pad->cd(1);
639 //gPad->SetBorderMode(0);
640 }
641
642 AppendPad(option);
643 //fGeomCam->AppendPad();
644
645 //
646 // Do not change gPad. The user should not see, that Draw
647 // changes gPad...
648 //
649// if (!hassame)
650// pad->cd();
651}
652
653// ------------------------------------------------------------------------
654//
655// Creates a TH1D which contains the projection of the contents of the
656// MHCamera onto the y-axis. The maximum and minimum are calculated
657// such that a slighly wider range than (GetMinimum(), GetMaximum()) is
658// displayed using THLimitsFinder::OptimizeLimits.
659//
660// If no name is given the newly allocated histogram is removed from
661// the current directory calling SetDirectory(0) in any other case
662// the newly created histogram is removed from the current directory
663// and added to gROOT such the gROOT->FindObject can find the histogram.
664//
665// If the standard name "_proj" is given "_proj" is appended to the name
666// of the MHCamera and the corresponding histogram is searched using
667// gROOT->FindObject and updated with the present projection.
668//
669// It is the responsibility of the user to make sure, that the newly
670// created histogram is freed correctly.
671//
672// Currently the new histogram is restrictred to 50 bins.
673// Maybe a optimal number can be calulated from the number of
674// bins on the x-axis of the MHCamera?
675//
676// The code was taken mainly from TH2::ProjectX such the interface
677// is more or less the same than to TH2-projections.
678//
679// If sector>=0 only entries with matching sector index are taken
680// into account.
681//
682TH1D *MHCamera::ProjectionS(const TArrayI &sector, const TArrayI &aidx, const char *name, const Int_t nbins) const
683{
684
685 // Create the projection histogram
686 TString pname(name);
687 if (pname=="_proj")
688 {
689 pname.Prepend(GetName());
690 if (sector.GetSize()>0)
691 {
692 pname += ";";
693 for (int i=0; i<sector.GetSize(); i++)
694 pname += sector[i];
695 }
696 if (aidx.GetSize()>0)
697 {
698 pname += ";";
699 for (int i=0; i<aidx.GetSize(); i++)
700 pname += aidx[i];
701 }
702 }
703
704 TH1D *h1=0;
705
706 //check if histogram with identical name exist
707 TObject *h1obj = gROOT->FindObject(pname);
708 if (h1obj && h1obj->InheritsFrom(TH1D::Class())) {
709 h1 = (TH1D*)h1obj;
710 h1->Reset();
711 }
712
713 if (!h1)
714 {
715 h1 = new TH1D;
716 h1->UseCurrentStyle();
717 h1->SetName(pname);
718 h1->SetTitle(GetTitle());
719 h1->SetDirectory(0);
720 h1->SetXTitle(GetYaxis()->GetTitle());
721 h1->SetYTitle("Counts");
722 //h1->Sumw2();
723 }
724
725 Double_t min = GetMinimumSectors(sector, aidx);
726 Double_t max = GetMaximumSectors(sector, aidx);
727
728 if (min==max && max>0)
729 min=0;
730 if (min==max && min<0)
731 max=0;
732
733 Int_t newbins=0;
734 THLimitsFinder::OptimizeLimits(nbins, newbins, min, max, kFALSE);
735
736 MBinning bins(nbins, min, max);
737 bins.Apply(*h1);
738
739 // Fill the projected histogram
740 for (Int_t idx=0; idx<fNcells-2; idx++)
741 if (IsUsed(idx) && MatchSector(idx, sector, aidx) && TMath::Finite(GetBinContent(idx+1)))
742 h1->Fill(GetBinContent(idx+1));
743
744 return h1;
745}
746
747// ------------------------------------------------------------------------
748//
749// Creates a TH1D which contains the projection of the contents of the
750// MHCamera onto the radius from the camera center.
751// The maximum and minimum are calculated
752// such that a slighly wider range than (GetMinimum(), GetMaximum()) is
753// displayed using THLimitsFinder::OptimizeLimits.
754//
755// If no name is given the newly allocated histogram is removed from
756// the current directory calling SetDirectory(0) in any other case
757// the newly created histogram is removed from the current directory
758// and added to gROOT such the gROOT->FindObject can find the histogram.
759//
760// If the standard name "_rad" is given "_rad" is appended to the name
761// of the MHCamera and the corresponding histogram is searched using
762// gROOT->FindObject and updated with the present projection.
763//
764// It is the responsibility of the user to make sure, that the newly
765// created histogram is freed correctly.
766//
767// Currently the new histogram is restrictred to 50 bins.
768// Maybe a optimal number can be calulated from the number of
769// bins on the x-axis of the MHCamera?
770//
771// The code was taken mainly from TH2::ProjectX such the interface
772// is more or less the same than to TH2-projections.
773//
774// If sector>=0 only entries with matching sector index are taken
775// into account.
776//
777TProfile *MHCamera::RadialProfileS(const TArrayI &sector, const TArrayI &aidx, const char *name, const Int_t nbins) const
778{
779 // Create the projection histogram
780 TString pname(name);
781 if (pname=="_rad")
782 {
783 pname.Prepend(GetName());
784 if (sector.GetSize()>0)
785 {
786 pname += ";";
787 for (int i=0; i<sector.GetSize(); i++)
788 pname += sector[i];
789 }
790 if (aidx.GetSize()>0)
791 {
792 pname += ";";
793 for (int i=0; i<aidx.GetSize(); i++)
794 pname += aidx[i];
795 }
796 }
797
798 TProfile *h1=0;
799
800 //check if histogram with identical name exist
801 TObject *h1obj = gROOT->FindObject(pname);
802 if (h1obj && h1obj->InheritsFrom(TProfile::Class())) {
803 h1 = (TProfile*)h1obj;
804 h1->Reset();
805 }
806
807 if (!h1)
808 {
809 h1 = new TProfile;
810 h1->UseCurrentStyle();
811 h1->SetName(pname);
812 h1->SetTitle(GetTitle());
813 h1->SetDirectory(0);
814 h1->SetXTitle("Radius from camera center [deg]");
815 h1->SetYTitle(GetYaxis()->GetTitle());
816 }
817
818 const Double_t m2d = fGeomCam->GetConvMm2Deg();
819
820 Double_t min = 0.;
821 Double_t max = fGeomCam->GetMaxRadius()*m2d;
822
823 Int_t newbins=0;
824
825 THLimitsFinder::OptimizeLimits(nbins, newbins, min, max, kFALSE);
826
827 MBinning bins(nbins, min, max);
828 bins.Apply(*h1);
829
830 // Fill the projected histogram
831 for (Int_t idx=0; idx<fNcells-2; idx++)
832 if (IsUsed(idx) && MatchSector(idx, sector, aidx) && TMath::Finite(GetBinContent(idx+1)))
833 h1->Fill(TMath::Hypot((*fGeomCam)[idx].GetX(),(*fGeomCam)[idx].GetY())*m2d,
834 GetBinContent(idx+1));
835 return h1;
836}
837
838
839// ------------------------------------------------------------------------
840//
841// Creates a TH1D which contains the projection of the contents of the
842// MHCamera onto the azimuth angle in the camera.
843//
844// If no name is given the newly allocated histogram is removed from
845// the current directory calling SetDirectory(0) in any other case
846// the newly created histogram is removed from the current directory
847// and added to gROOT such the gROOT->FindObject can find the histogram.
848//
849// If the standard name "_az" is given "_az" is appended to the name
850// of the MHCamera and the corresponding histogram is searched using
851// gROOT->FindObject and updated with the present projection.
852//
853// It is the responsibility of the user to make sure, that the newly
854// created histogram is freed correctly.
855//
856// Currently the new histogram is restrictred to 60 bins.
857// Maybe a optimal number can be calulated from the number of
858// bins on the x-axis of the MHCamera?
859//
860// The code was taken mainly from TH2::ProjectX such the interface
861// is more or less the same than to TH2-projections.
862//
863TProfile *MHCamera::AzimuthProfileA(const TArrayI &aidx, const char *name, const Int_t nbins) const
864{
865 // Create the projection histogram
866 TString pname(name);
867 if (pname=="_az")
868 {
869 pname.Prepend(GetName());
870 if (aidx.GetSize()>0)
871 {
872 pname += ";";
873 for (int i=0; i<aidx.GetSize(); i++)
874 pname += aidx[i];
875 }
876 }
877
878 TProfile *h1=0;
879
880 //check if histogram with identical name exist
881 TObject *h1obj = gROOT->FindObject(pname);
882 if (h1obj && h1obj->InheritsFrom(TProfile::Class())) {
883 h1 = (TProfile*)h1obj;
884 h1->Reset();
885 }
886
887 if (!h1)
888 {
889
890 h1 = new TProfile;
891 h1->UseCurrentStyle();
892 h1->SetName(pname);
893 h1->SetTitle(GetTitle());
894 h1->SetDirectory(0);
895 h1->SetXTitle("Azimuth in camera [deg]");
896 h1->SetYTitle(GetYaxis()->GetTitle());
897 }
898
899 //Double_t min = 0;
900 //Double_t max = 360;
901
902 //Int_t newbins=0;
903 //THLimitsFinder::OptimizeLimits(nbins, newbins, min, max, kFALSE);
904
905 MBinning bins(nbins, 0, 360);
906 bins.Apply(*h1);
907
908 // Fill the projected histogram
909 for (Int_t idx=0; idx<fNcells-2; idx++)
910 {
911 if (IsUsed(idx) && MatchSector(idx, TArrayI(), aidx) && TMath::Finite(GetPixContent(idx)))
912 h1->Fill(TMath::ATan2((*fGeomCam)[idx].GetY(),(*fGeomCam)[idx].GetX())*TMath::RadToDeg()+180,
913 GetPixContent(idx));
914
915 }
916
917 return h1;
918}
919
920// ------------------------------------------------------------------------
921//
922// Updates the pixel colors and paints the pixels
923//
924void MHCamera::Update(Bool_t islog, Bool_t isbox, Bool_t iscol, Bool_t issame)
925{
926 Double_t min = GetMinimum(kFALSE);
927 Double_t max = GetMaximum(kFALSE);
928 if (min==FLT_MAX)
929 {
930 min = 0;
931 max = 1;
932 }
933
934 if (min==max)
935 max += 1;
936
937 if (!issame)
938 UpdateLegend(min, max, islog);
939
940 // Try to estimate the units of the current display. This is only
941 // necessary for 'same' option and allows distorted images of the camera!
942 const Float_t maxr = (1-fGeomCam->GetConvMm2Deg())*fGeomCam->GetMaxRadius()/2;
943 const Float_t conv = !issame ||
944 gPad->GetX1()<-maxr || gPad->GetY1()<-maxr ||
945 gPad->GetX2()> maxr || gPad->GetY2()>maxr ? 1 : fGeomCam->GetConvMm2Deg();
946
947 TAttLine line(kBlack, kSolid, 1);
948 TAttFill fill;
949 for (Int_t i=0; i<fNcells-2; i++)
950 {
951 fill.SetFillStyle(issame || (IsTransparent() && !IsUsed(i)) ? 0 : 1001);
952
953 if (!issame)
954 {
955 const Bool_t isnan = !TMath::Finite(fArray[i+1]);
956 if (!IsUsed(i) || !iscol || isnan)
957 {
958 fill.SetFillColor(10);
959
960 if (isnan)
961 gLog << warn << "MHCamera::Update: " << GetName() << " <" << GetTitle() << "> - Pixel Index #" << i << " contents is not finite..." << endl;
962 }
963 else
964 fill.SetFillColor(GetColor(GetBinContent(i+1), min, max, islog));
965 }
966
967 const MGeom &pix = (*fGeomCam)[i];
968
969 Double_t scale = 1;//conv/(fAbberation+1);
970
971 if (!isbox && !IsUsed(i) && TestBit(kNoUnused))
972 continue;
973
974 if (isbox && (!IsUsed(i) || !TMath::Finite(fArray[i+1])))
975 continue;
976
977 if (isbox)
978 {
979 scale = (GetBinContent(i+1)-min)/(max-min);
980 if (scale>1)
981 scale=1;
982 }
983
984 pix.PaintPrimitive(line, fill, conv, scale/(fAbberation+1));
985 }
986}
987
988// ------------------------------------------------------------------------
989//
990// Print minimum and maximum
991//
992void MHCamera::Print(Option_t *) const
993{
994 gLog << all << "Minimum: " << GetMinimum();
995 if (fMinimum==-1111)
996 gLog << " <autoscaled>";
997 gLog << endl;
998 gLog << "Maximum: " << GetMaximum();
999 if (fMaximum==-1111)
1000 gLog << " <autoscaled>";
1001 gLog << endl;
1002}
1003
1004// ------------------------------------------------------------------------
1005//
1006// Paint the y-axis title
1007//
1008void MHCamera::PaintAxisTitle()
1009{
1010 const Float_t range = fGeomCam->GetMaxRadius()*1.05;
1011 const Float_t w = (1 + 1.5/sqrt((float)(fNcells-2)))*range;
1012
1013 TLatex *ptitle = new TLatex(w, -.90*range, GetYaxis()->GetTitle());
1014
1015 ptitle->SetTextSize(0.05);
1016 ptitle->SetTextAlign(21);
1017
1018 // box with the histogram title
1019 ptitle->SetTextColor(gStyle->GetTitleTextColor());
1020#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,01)
1021 ptitle->SetTextFont(gStyle->GetTitleFont(""));
1022#endif
1023 ptitle->Paint();
1024}
1025
1026// ------------------------------------------------------------------------
1027//
1028// Paints the camera.
1029//
1030void MHCamera::Paint(Option_t *o)
1031{
1032 if (fNcells<=1)
1033 return;
1034
1035 TString opt(o);
1036 opt.ToLower();
1037
1038 if (opt.Contains("hist"))
1039 {
1040 opt.ReplaceAll("hist", "");
1041 opt.ReplaceAll("box", "");
1042 opt.ReplaceAll("pixelindex", "");
1043 opt.ReplaceAll("sectorindex", "");
1044 opt.ReplaceAll("abscontent", "");
1045 opt.ReplaceAll("content", "");
1046 opt.ReplaceAll("integer", "");
1047 opt.ReplaceAll("float", "");
1048 opt.ReplaceAll("proj", "");
1049 opt.ReplaceAll("pal0", "");
1050 opt.ReplaceAll("pal1", "");
1051 opt.ReplaceAll("pal2", "");
1052 opt.ReplaceAll("nopal", "");
1053 TH1D::Paint(opt);
1054 return;
1055 }
1056
1057 if (opt.Contains("proj"))
1058 {
1059 opt.ReplaceAll("proj", "");
1060 Projection(GetName())->Paint(opt);
1061 return;
1062 }
1063
1064 const Bool_t hassame = opt.Contains("same");
1065 const Bool_t hasbox = opt.Contains("box");
1066 const Bool_t hascol = hasbox ? !opt.Contains("nocol") : kTRUE;
1067
1068 if (!hassame)
1069 gPad->Clear();
1070
1071 if (!fGeomCam)
1072 return;
1073
1074 if (!hassame)
1075 {
1076 // Maintain aspect ratio
1077 const Float_t r = fGeomCam->GetMaxRadius();
1078
1079 MH::SetPadRange(-r*1.02, -r*1.02, TestBit(kNoLegend) ? r : 1.15*r, r*1.2);
1080
1081 if (GetPainter())
1082 {
1083 // Paint statistics
1084 if (!TestBit(TH1::kNoStats))
1085 fPainter->PaintStat(gStyle->GetOptStat(), NULL);
1086
1087 // Paint primitives (pixels, color legend, photons, ...)
1088 if (fPainter->InheritsFrom(THistPainter::Class()))
1089 {
1090 static_cast<THistPainter*>(fPainter)->MakeChopt("");
1091 static_cast<THistPainter*>(fPainter)->PaintTitle();
1092 }
1093 }
1094 }
1095
1096 const Bool_t pal1 = opt.Contains("pal1");
1097 const Bool_t pal2 = opt.Contains("pal2");
1098 const Bool_t nopal = opt.Contains("nopal");
1099
1100 if (!pal1 && !pal2 && !nopal)
1101 SetPrettyPalette();
1102
1103 if (pal1)
1104 SetDeepBlueSeaPalette();
1105
1106 if (pal2)
1107 SetInvDeepBlueSeaPalette();
1108
1109 // Update Contents of the pixels and paint legend
1110 Update(gPad->GetLogy(), hasbox, hascol, hassame);
1111
1112 if (!hassame)
1113 PaintAxisTitle();
1114
1115 if (opt.Contains("pixelindex"))
1116 {
1117 PaintIndices(0);
1118 return;
1119 }
1120 if (opt.Contains("sectorindex"))
1121 {
1122 PaintIndices(1);
1123 return;
1124 }
1125 if (opt.Contains("abscontent"))
1126 {
1127 PaintIndices(3);
1128 return;
1129 }
1130 if (opt.Contains("content"))
1131 {
1132 PaintIndices(2);
1133 return;
1134 }
1135 if (opt.Contains("pixelentries"))
1136 {
1137 PaintIndices(4);
1138 return;
1139 }
1140 if (opt.Contains("text"))
1141 {
1142 PaintIndices(5);
1143 return;
1144 }
1145 if (opt.Contains("integer"))
1146 {
1147 PaintIndices(6);
1148 return;
1149 }
1150 if (opt.Contains("float"))
1151 {
1152 PaintIndices(7);
1153 return;
1154 }
1155}
1156
1157void MHCamera::SetDrawOption(Option_t *option)
1158{
1159 // This is a workaround. For some reason MHCamera is
1160 // stored in a TObjLink instead of a TObjOptLink
1161 if (!option || !gPad)
1162 return;
1163
1164 TListIter next(gPad->GetListOfPrimitives());
1165 delete gPad->FindObject("Tframe");
1166 TObject *obj;
1167 while ((obj = next()))
1168 if (obj == this && (TString)next.GetOption()!=(TString)option)
1169 {
1170 gPad->GetListOfPrimitives()->Remove(this);
1171 gPad->GetListOfPrimitives()->AddFirst(this, option);
1172 return;
1173 }
1174}
1175
1176// ------------------------------------------------------------------------
1177//
1178// With this function you can change the color palette. For more
1179// information see TStyle::SetPalette. Only palettes with 50 colors
1180// are allowed.
1181// In addition you can use SetPalette(52, 0) to create an inverse
1182// deep blue sea palette
1183//
1184void MHCamera::SetPalette(Int_t ncolors, Int_t *colors)
1185{
1186 //
1187 // If not enough colors are specified skip this.
1188 //
1189 if (ncolors>1 && ncolors<50)
1190 {
1191 gLog << err << "MHCamera::SetPalette: Only default palettes with 50 colors are allowed... ignored." << endl;
1192 return;
1193 }
1194
1195 //
1196 // If ncolors==52 create a reversed deep blue sea palette
1197 //
1198 if (ncolors==52)
1199 {
1200 gStyle->SetPalette(51, NULL);
1201
1202 const Int_t n = GetContour()==0?50:GetContour();
1203
1204 TArrayI c(n);
1205 for (int i=0; i<n; i++)
1206 c[n-i-1] = gStyle->GetColorPalette(i);
1207 gStyle->SetPalette(n, c.GetArray());
1208 }
1209 else
1210 gStyle->SetPalette(ncolors, colors);
1211}
1212
1213
1214// ------------------------------------------------------------------------
1215//
1216// Changes the palette of the displayed camera histogram.
1217//
1218// Change to the right pad first - otherwise GetDrawOption() might fail.
1219//
1220void MHCamera::SetPrettyPalette()
1221{
1222 TString opt(GetDrawOption());
1223
1224 if (!opt.Contains("hist", TString::kIgnoreCase))
1225 SetPalette(1, 0);
1226
1227 opt.ReplaceAll("pal1", "");
1228 opt.ReplaceAll("pal2", "");
1229
1230 SetDrawOption(opt);
1231}
1232
1233// ------------------------------------------------------------------------
1234//
1235// Changes the palette of the displayed camera histogram.
1236//
1237// Change to the right pad first - otherwise GetDrawOption() might fail.
1238//
1239void MHCamera::SetDeepBlueSeaPalette()
1240{
1241 TString opt(GetDrawOption());
1242
1243 if (!opt.Contains("hist", TString::kIgnoreCase))
1244 SetPalette(51, 0);
1245
1246 opt.ReplaceAll("pal1", "");
1247 opt.ReplaceAll("pal2", "");
1248 opt += "pal1";
1249
1250 SetDrawOption(opt);
1251}
1252
1253// ------------------------------------------------------------------------
1254//
1255// Changes the palette of the displayed camera histogram.
1256//
1257// Change to the right pad first - otherwise GetDrawOption() might fail.
1258//
1259void MHCamera::SetInvDeepBlueSeaPalette()
1260{
1261 TString opt(GetDrawOption());
1262
1263 if (!opt.Contains("hist", TString::kIgnoreCase))
1264 SetPalette(52, 0);
1265
1266 opt.ReplaceAll("pal1", "");
1267 opt.ReplaceAll("pal2", "");
1268 opt += "pal2";
1269
1270 SetDrawOption(opt);
1271}
1272
1273// ------------------------------------------------------------------------
1274//
1275// Paint indices (as text) inside the pixels. Depending of the type-
1276// argument we paint:
1277// 0: pixel number
1278// 1: sector number
1279// 2: content
1280// 5: Assume GetBinContent is a char
1281//
1282void MHCamera::PaintIndices(Int_t type)
1283{
1284 if (fNcells<=1)
1285 return;
1286
1287 const Double_t min = GetMinimum();
1288 const Double_t max = GetMaximum();
1289
1290 if (type==2 && max==min)
1291 return;
1292
1293 TText txt;
1294 if (type!=5)
1295 txt.SetTextFont(122);
1296 txt.SetTextAlign(22); // centered/centered
1297
1298 for (Int_t i=0; i<fNcells-2; i++)
1299 {
1300 const MGeom &h = (*fGeomCam)[i];
1301
1302 TString num;
1303 switch (type)
1304 {
1305 case 0: num += i; break;
1306 case 1: num += h.GetSector(); break;
1307 case 2: num += TMath::Nint((fArray[i+1]-min)/(max-min)); break;
1308 case 3: num += TMath::Nint(fArray[i+1]); break;
1309 case 4: num += fBinEntries[i+1]; break;
1310 case 5: num = (char)TMath::Nint(GetBinContent(i+1)); break;
1311 case 6: num += TMath::Nint(GetBinContent(i+1)); break;
1312 case 7: num += GetBinContent(i+1); break;
1313 }
1314
1315 // FIXME: Should depend on the color of the pixel...
1316 //(GetColor(GetBinContent(i+1), min, max, 0));
1317 txt.SetTextColor(kRed);
1318 txt.SetTextSize(0.3*h.GetT()/fGeomCam->GetMaxRadius()/1.05);
1319 txt.PaintText(h.GetX(), h.GetY(), num);
1320 }
1321}
1322
1323// ------------------------------------------------------------------------
1324//
1325// Call this function to add a MCamEvent on top of the present contents.
1326//
1327void MHCamera::AddMeanShift(const MCamEvent &event, Int_t type, Stat_t w)
1328{
1329 if (fNcells<=1 || IsFreezed())
1330 return;
1331
1332 const Double_t mean = event.GetCameraMean(*fGeomCam, type);
1333
1334 // FIXME: Security check missing!
1335 for (Int_t idx=0; idx<fNcells-2; idx++)
1336 {
1337 Double_t val=0;
1338 if (event.GetPixelContent(val, idx, *fGeomCam, type)/* && !IsUsed(idx)*/)
1339 {
1340 SetUsed(idx);
1341 Fill(idx, (val-mean)*w); // FIXME: Slow!
1342 }
1343 }
1344 fEntries++;
1345}
1346
1347// ------------------------------------------------------------------------
1348//
1349// Call this function to add a MCamEvent on top of the present contents.
1350//
1351void MHCamera::AddMedianShift(const MCamEvent &event, Int_t type, Stat_t w)
1352{
1353 if (fNcells<=1 || IsFreezed())
1354 return;
1355
1356 const Double_t median = event.GetCameraMedian(*fGeomCam, type);
1357
1358 // FIXME: Security check missing!
1359 for (Int_t idx=0; idx<fNcells-2; idx++)
1360 {
1361 Double_t val=0;
1362 if (event.GetPixelContent(val, idx, *fGeomCam, type)/* && !IsUsed(idx)*/)
1363 {
1364 SetUsed(idx);
1365 Fill(idx, (val-median)*w); // FIXME: Slow!
1366 }
1367 }
1368 fEntries++;
1369}
1370
1371// ------------------------------------------------------------------------
1372//
1373// Call this function to add a MCamEvent on top of the present contents.
1374//
1375void MHCamera::AddCamContent(const MCamEvent &event, Int_t type, Stat_t w)
1376{
1377 if (fNcells<=1 || IsFreezed())
1378 return;
1379
1380 // FIXME: Security check missing!
1381 for (Int_t idx=0; idx<fNcells-2; idx++)
1382 {
1383 Double_t val=0;
1384 if (event.GetPixelContent(val, idx, *fGeomCam, type)/* && !IsUsed(idx)*/)
1385 {
1386 SetUsed(idx);
1387 Fill(idx, val*w); // FIXME: Slow!
1388 }
1389 }
1390 fEntries++;
1391}
1392
1393// ------------------------------------------------------------------------
1394//
1395void MHCamera::AddCamDifference(const MCamEvent &event1, const MCamEvent &event2, Int_t type, Stat_t weight)
1396{
1397 if (fNcells<=1 || IsFreezed())
1398 return;
1399
1400 // FIXME: Security check missing!
1401 for (Int_t idx=0; idx<fNcells-2; idx++)
1402 {
1403 Double_t val1=0;
1404 if (!event1.GetPixelContent(val1, idx, *fGeomCam, type))
1405 continue;
1406
1407 Double_t val2=0;
1408 if (!event2.GetPixelContent(val2, idx, *fGeomCam, type))
1409 continue;
1410
1411 SetUsed(idx);
1412 Fill(idx, (val1-val2)*weight); // FIXME: Slow!
1413 }
1414 fEntries++;
1415}
1416
1417// ------------------------------------------------------------------------
1418//
1419void MHCamera::AddCamRatio(const MCamEvent &event1, const MCamEvent &event2, Int_t type, Stat_t weight)
1420{
1421 if (fNcells<=1 || IsFreezed())
1422 return;
1423
1424 // FIXME: Security check missing!
1425 for (Int_t idx=0; idx<fNcells-2; idx++)
1426 {
1427 Double_t val1=0;
1428 if (!event1.GetPixelContent(val1, idx, *fGeomCam, type))
1429 continue;
1430
1431 Double_t val2=0;
1432 if (!event2.GetPixelContent(val2, idx, *fGeomCam, type))
1433 continue;
1434
1435 if (val2==0)
1436 continue;
1437
1438 SetUsed(idx);
1439 Fill(idx, weight*val1/val2); // FIXME: Slow!
1440 }
1441 fEntries++;
1442}
1443
1444// ------------------------------------------------------------------------
1445//
1446// Call this function to add a MCamEvent on top of the present contents.
1447//
1448void MHCamera::SetCamError(const MCamEvent &evt, Int_t type)
1449{
1450
1451 if (fNcells<=1 || IsFreezed())
1452 return;
1453
1454 // FIXME: Security check missing!
1455 for (Int_t idx=0; idx<fNcells-2; idx++)
1456 {
1457 Double_t val=0;
1458 if (evt.GetPixelContent(val, idx, *fGeomCam, type)/* && !IsUsed(idx)*/)
1459 SetUsed(idx);
1460
1461 SetBinError(idx+1, val); // FIXME: Slow!
1462 }
1463}
1464
1465Stat_t MHCamera::GetBinContent(Int_t bin) const
1466{
1467 if (fBuffer) ((TH1D*)this)->BufferEmpty();
1468 if (bin < 0) bin = 0;
1469 if (bin >= fNcells) bin = fNcells-1;
1470 if (!fArray) return 0;
1471
1472 if (!TestBit(kProfile))
1473 return Stat_t (fArray[bin]);
1474
1475 if (fBinEntries.fArray[bin] == 0) return 0;
1476 return fArray[bin]/fBinEntries.fArray[bin];
1477}
1478
1479// ------------------------------------------------------------------------
1480//
1481// In the case the kProfile flag is set the spread of the bin is returned.
1482// If you want to have the mean error instead set the kErrorMean bit via
1483// SetBit(kErrorMean) first.
1484//
1485Stat_t MHCamera::GetBinError(Int_t bin) const
1486{
1487 if (!TestBit(kProfile))
1488 return TH1D::GetBinError(bin);
1489
1490 const UInt_t n = (UInt_t)fBinEntries[bin];
1491
1492 if (n==0)
1493 return 0;
1494
1495 const Double_t sqr = fSumw2.fArray[bin] / n;
1496 const Double_t val = fArray[bin] / n;
1497
1498 const Double_t spread = sqr>val*val ? TMath::Sqrt(sqr - val*val) : 0;
1499
1500 return TestBit(kErrorMean) ? spread/TMath::Sqrt(n) : spread;
1501
1502 /*
1503 Double_t rc = 0;
1504 if (TestBit(kSqrtVariance) && GetEntries()>0) // error on the mean
1505 {
1506 const Double_t error = fSumw2.fArray[bin]/GetEntries();
1507 const Double_t val = fArray[bin]/GetEntries();
1508 rc = val*val>error ? 0 : TMath::Sqrt(error - val*val);
1509 }
1510 else
1511 rc = TH1D::GetBinError(bin);
1512
1513 return Profile(rc);*/
1514}
1515
1516// ------------------------------------------------------------------------
1517//
1518// Call this function to add a MHCamera on top of the present contents.
1519// Type:
1520// 0) bin content
1521// 1) errors
1522// 2) rel. errors
1523//
1524void MHCamera::AddCamContent(const MHCamera &d, Int_t type)
1525{
1526 if (fNcells!=d.fNcells || IsFreezed())
1527 return;
1528
1529 // FIXME: Security check missing!
1530 for (Int_t idx=0; idx<fNcells-2; idx++)
1531 if (d.IsUsed(idx))
1532 SetUsed(idx);
1533
1534 switch (type)
1535 {
1536 case 1:
1537 // Under-/Overflow bins not handled!
1538 for (Int_t idx=0; idx<fNcells-2; idx++)
1539 if (d.IsUsed(idx))
1540 Fill(idx, d.GetBinError(idx+1));
1541 fEntries++;
1542 break;
1543 case 2:
1544 // Under-/Overflow bins not handled!
1545 for (Int_t idx=0; idx<fNcells-2; idx++)
1546 if (d.GetBinContent(idx+1)!=0 && d.IsUsed(idx))
1547 Fill(idx, TMath::Abs(d.GetBinError(idx+1)/d.GetBinContent(idx+1)));
1548 fEntries++;
1549 break;
1550 default:
1551 if (TestBit(kProfile)!=d.TestBit(kProfile))
1552 gLog << warn << "WARNING - You have tried to call AddCamContent for two different kind of histograms (kProfile set or not)." << endl;
1553
1554 // environment
1555 fEntries += d.fEntries;
1556 fTsumw += d.fTsumw;
1557 fTsumw2 += d.fTsumw2;
1558 fTsumwx += d.fTsumwx;
1559 fTsumwx2 += d.fTsumwx2;
1560 // Bin contents
1561 for (Int_t idx=1; idx<fNcells-1; idx++)
1562 {
1563 if (!d.IsUsed(idx-1))
1564 continue;
1565
1566 fArray[idx] += d.fArray[idx];
1567 fBinEntries[idx] += d.fBinEntries[idx];
1568 fSumw2.fArray[idx] += d.fSumw2.fArray[idx];
1569 }
1570 // Underflow bin
1571 fArray[0] += d.fArray[0];
1572 fBinEntries[0] += d.fBinEntries[0];
1573 fSumw2.fArray[0] += d.fSumw2.fArray[0];
1574 // Overflow bin
1575 fArray[fNcells-1] += d.fArray[fNcells-1];
1576 fBinEntries[fNcells-1] += d.fBinEntries[fNcells-1];
1577 fSumw2.fArray[fNcells-1] += d.fSumw2.fArray[fNcells-1];
1578 break;
1579/* default:
1580 if (TestBit(kProfile)!=d.TestBit(kProfile))
1581 gLog << warn << "WARNING - You have tried to call AddCamContent for two different kind of histograms (kProfile set or not)." << endl;
1582
1583 for (Int_t idx=0; idx<fNcells-2; idx++)
1584 Fill(idx, d.GetBinContent(idx+1));
1585 break;*/
1586 }
1587 fEntries++;
1588}
1589
1590// ------------------------------------------------------------------------
1591//
1592// Call this function to add a TArrayD on top of the present contents.
1593//
1594void MHCamera::AddCamContent(const TArrayD &event, const TArrayC *used)
1595{
1596 if (event.GetSize()!=fNcells-2 || IsFreezed())
1597 return;
1598
1599 if (used && used->GetSize()!=fNcells-2)
1600 return;
1601
1602 for (Int_t idx=0; idx<fNcells-2; idx++)
1603 {
1604 Fill(idx, event[idx]); // FIXME: Slow!
1605
1606 if (!used || (*used)[idx])
1607 SetUsed(idx);
1608 }
1609 fEntries++;
1610}
1611
1612// ------------------------------------------------------------------------
1613//
1614// Call this function to add a MArrayD on top of the present contents.
1615//
1616void MHCamera::AddCamContent(const MArrayD &event, const TArrayC *used)
1617{
1618 if (event.GetSize()!=(UInt_t)(fNcells-2) || IsFreezed())
1619 return;
1620
1621 if (used && used->GetSize()!=fNcells-2)
1622 return;
1623
1624 for (Int_t idx=0; idx<fNcells-2; idx++)
1625 {
1626 Fill(idx, event[idx]); // FIXME: Slow!
1627
1628 if (!used || (*used)[idx])
1629 SetUsed(idx);
1630 }
1631 fEntries++;
1632}
1633
1634// ------------------------------------------------------------------------
1635//
1636// Call this function to add a MCamEvent on top of the present contents.
1637// 1 is added to each pixel if the contents of MCamEvent>threshold
1638// (in case isabove is set to kTRUE == default)
1639// 1 is added to each pixel if the contents of MCamEvent<threshold
1640// (in case isabove is set to kFALSE)
1641//
1642// in unused pixel is not counted if it didn't fullfill the condition.
1643//
1644void MHCamera::CntCamContent(const MCamEvent &event, Double_t threshold, Int_t type, Bool_t isabove)
1645{
1646 if (fNcells<=1 || IsFreezed())
1647 return;
1648
1649 // FIXME: Security check missing!
1650 for (Int_t idx=0; idx<fNcells-2; idx++)
1651 {
1652 Double_t val=threshold;
1653 const Bool_t rc = event.GetPixelContent(val, idx, *fGeomCam, type);
1654 if (rc)
1655 SetUsed(idx);
1656
1657 const Bool_t cond =
1658 ( isabove && val>threshold) ||
1659 (!isabove && val<threshold);
1660
1661 Fill(idx, rc && cond ? 1 : 0);
1662 }
1663 fEntries++;
1664}
1665
1666// ------------------------------------------------------------------------
1667//
1668// Call this function to add a MCamEvent on top of the present contents.
1669// - the contents of the pixels in event are added to each pixel
1670// if the pixel of thresevt<threshold (in case isabove is set
1671// to kTRUE == default)
1672// - the contents of the pixels in event are added to each pixel
1673// if the pixel of thresevt<threshold (in case isabove is set
1674// to kFALSE)
1675//
1676// in unused pixel is not counted if it didn't fullfill the condition.
1677//
1678void MHCamera::CntCamContent(const MCamEvent &event, Int_t type1, const MCamEvent &thresevt, Int_t type2, Double_t threshold, Bool_t isabove)
1679{
1680 if (fNcells<=1 || IsFreezed())
1681 return;
1682
1683 // FIXME: Security check missing!
1684 for (Int_t idx=0; idx<fNcells-2; idx++)
1685 {
1686 Double_t th=0;
1687 if (!thresevt.GetPixelContent(th, idx, *fGeomCam, type2))
1688 continue;
1689
1690 if ((isabove && th>threshold) || (!isabove && th<threshold))
1691 continue;
1692
1693 Double_t val=th;
1694 if (event.GetPixelContent(val, idx, *fGeomCam, type1))
1695 {
1696 SetUsed(idx);
1697 Fill(idx, val);
1698 }
1699 }
1700 fEntries++;
1701}
1702
1703// ------------------------------------------------------------------------
1704//
1705// Call this function to add a MCamEvent on top of the present contents.
1706// 1 is added to each pixel if the contents of MCamEvent>threshold
1707// (in case isabove is set to kTRUE == default)
1708// 1 is added to each pixel if the contents of MCamEvent<threshold
1709// (in case isabove is set to kFALSE)
1710//
1711// in unused pixel is not counted if it didn't fullfill the condition.
1712//
1713void MHCamera::CntCamContent(const MCamEvent &event, TArrayD threshold, Int_t type, Bool_t isabove)
1714{
1715 if (fNcells<=1 || IsFreezed())
1716 return;
1717
1718 // FIXME: Security check missing!
1719 for (Int_t idx=0; idx<fNcells-2; idx++)
1720 {
1721 Double_t val=threshold[idx];
1722 if (event.GetPixelContent(val, idx, *fGeomCam, type)/* && !IsUsed(idx)*/)
1723 {
1724 SetUsed(idx);
1725
1726 if (val>threshold[idx] && isabove)
1727 Fill(idx);
1728 if (val<threshold[idx] && !isabove)
1729 Fill(idx);
1730 }
1731 }
1732 fEntries++;
1733}
1734
1735// ------------------------------------------------------------------------
1736//
1737// Call this function to add a TArrayD on top of the present contents.
1738// 1 is added to each pixel if the contents of MCamEvent>threshold
1739//
1740void MHCamera::CntCamContent(const TArrayD &event, Double_t threshold, Bool_t ispos)
1741{
1742 if (event.GetSize()!=fNcells-2 || IsFreezed())
1743 return;
1744
1745 for (Int_t idx=0; idx<fNcells-2; idx++)
1746 {
1747 if (event[idx]>threshold)
1748 Fill(idx);
1749
1750 if (!ispos || fArray[idx+1]>0)
1751 SetUsed(idx);
1752 }
1753 fEntries++;
1754}
1755
1756// ------------------------------------------------------------------------
1757//
1758// Call this function to add a MCamEvent on top of the present contents.
1759// 1 is added to each pixel if the contents of MCamEvent>threshold
1760// (in case isabove is set to kTRUE == default)
1761// 1 is added to each pixel if the contents of MCamEvent<threshold
1762// (in case isabove is set to kFALSE)
1763//
1764// in unused pixel is not counted if it didn't fullfill the condition.
1765//
1766void MHCamera::SetMaxCamContent(const MCamEvent &event, Int_t type)
1767{
1768 if (fNcells<=1 || IsFreezed())
1769 return;
1770
1771 // FIXME: Security check missing!
1772 for (Int_t idx=0; idx<fNcells-2; idx++)
1773 {
1774 Double_t val=0;
1775 const Bool_t rc = event.GetPixelContent(val, idx, *fGeomCam, type);
1776 if (!rc)
1777 continue;
1778
1779 if (!IsUsed(idx))
1780 {
1781 fArray[idx+1] = val;
1782 SetUsed(idx);
1783 fBinEntries.fArray[idx+1]=1;
1784 }
1785 else
1786 if (val>fArray[idx+1])
1787 fArray[idx+1] = val;
1788 }
1789 fEntries++;
1790}
1791
1792// ------------------------------------------------------------------------
1793//
1794// Call this function to add a MCamEvent on top of the present contents.
1795// 1 is added to each pixel if the contents of MCamEvent>threshold
1796// (in case isabove is set to kTRUE == default)
1797// 1 is added to each pixel if the contents of MCamEvent<threshold
1798// (in case isabove is set to kFALSE)
1799//
1800// in unused pixel is not counted if it didn't fullfill the condition.
1801//
1802void MHCamera::SetMinCamContent(const MCamEvent &event, Int_t type)
1803{
1804 if (fNcells<=1 || IsFreezed())
1805 return;
1806
1807 // FIXME: Security check missing!
1808 for (Int_t idx=0; idx<fNcells-2; idx++)
1809 {
1810 Double_t val=0;
1811 const Bool_t rc = event.GetPixelContent(val, idx, *fGeomCam, type);
1812 if (!rc)
1813 continue;
1814
1815 if (!IsUsed(idx))
1816 {
1817 fArray[idx+1] = val;
1818 SetUsed(idx);
1819 fBinEntries.fArray[idx+1]=1;
1820 }
1821 else
1822 if (val<fArray[idx+1])
1823 fArray[idx+1] = val;
1824 }
1825 fEntries++;
1826}
1827
1828// ------------------------------------------------------------------------
1829//
1830// Fill the pixels with random contents.
1831//
1832void MHCamera::FillRandom()
1833{
1834 if (fNcells<=1 || IsFreezed())
1835 return;
1836
1837 Reset();
1838
1839 // FIXME: Security check missing!
1840 for (Int_t idx=0; idx<fNcells-2; idx++)
1841 {
1842 Fill(idx, gRandom->Uniform()*fGeomCam->GetPixRatio(idx));
1843 SetUsed(idx);
1844 }
1845 fEntries=1;
1846}
1847
1848
1849// ------------------------------------------------------------------------
1850//
1851// The array must be in increasing order, eg: 2.5, 3.7, 4.9
1852// The values in each bin are replaced by the interval in which the value
1853// fits. In the example we have four intervals
1854// (<2.5, 2.5-3.7, 3.7-4.9, >4.9). Maximum and minimum are set
1855// accordingly.
1856//
1857void MHCamera::SetLevels(const TArrayF &arr)
1858{
1859 if (fNcells<=1)
1860 return;
1861
1862 for (Int_t i=0; i<fNcells-2; i++)
1863 {
1864 if (!IsUsed(i))
1865 continue;
1866
1867 Int_t j = arr.GetSize();
1868 while (j && fArray[i+1]<arr[j-1])
1869 j--;
1870
1871 fArray[i+1] = j;
1872 }
1873 SetMaximum(arr.GetSize());
1874 SetMinimum(0);
1875}
1876
1877// ------------------------------------------------------------------------
1878//
1879// Reset the all pixel colors to a default value
1880//
1881void MHCamera::Reset(Option_t *opt)
1882{
1883 if (fNcells<=1 || IsFreezed())
1884 return;
1885
1886 TH1::Reset(opt);
1887
1888 fUsed.Reset();
1889 fBinEntries.Reset();
1890
1891 for (Int_t i=0; i<fNcells; i++)
1892 fArray[i] = 0;
1893}
1894
1895// ------------------------------------------------------------------------
1896//
1897// Here we calculate the color index for the current value.
1898// The color index is defined with the class TStyle and the
1899// Color palette inside. We use the command gStyle->SetPalette(1,0)
1900// for the display. So we have to convert the value "wert" into
1901// a color index that fits the color palette.
1902// The range of the color palette is defined by the values fMinPhe
1903// and fMaxRange. Between this values we have 50 color index, starting
1904// with 0 up to 49.
1905//
1906Int_t MHCamera::GetColor(Float_t val, Float_t min, Float_t max, Bool_t islog)
1907{
1908 //
1909 // first treat the over- and under-flows
1910 //
1911 const Int_t ncol = GetContour()==0?50:GetContour();
1912
1913 const Int_t maxcolidx = ncol-1;
1914
1915 if (!TMath::Finite(val)) // FIXME: gLog!
1916 return maxcolidx/2;
1917
1918 if (val >= max)
1919 return gStyle->GetColorPalette(maxcolidx);
1920
1921 if (val <= min)
1922 return gStyle->GetColorPalette(0);
1923
1924 //
1925 // calculate the color index
1926 //
1927 Float_t ratio;
1928 if (islog && min>0)
1929 ratio = log10(val/min) / log10(max/min);
1930 else
1931 ratio = (val-min) / (max-min);
1932
1933 const Int_t colidx = TMath::FloorNint(ratio*ncol);
1934 return gStyle->GetColorPalette(colidx);
1935}
1936
1937TPaveStats *MHCamera::GetStatisticBox()
1938{
1939 TObject *obj = 0;
1940
1941 TIter Next(fFunctions);
1942 while ((obj = Next()))
1943 if (obj->InheritsFrom(TPaveStats::Class()))
1944 return static_cast<TPaveStats*>(obj);
1945
1946 return NULL;
1947}
1948
1949// ------------------------------------------------------------------------
1950//
1951// Change the text on the legend according to the range of the Display
1952//
1953void MHCamera::UpdateLegend(Float_t min, Float_t max, Bool_t islog)
1954{
1955 const Float_t range = fGeomCam->GetMaxRadius();
1956
1957 if (!TestBit(kNoScale))
1958 {
1959 TArrow arr;
1960 arr.PaintArrow(-range*.99, -range*.99, -range*.60, -range*.99, 0.010);
1961 arr.PaintArrow(-range*.99, -range*.99, -range*.99, -range*.60, 0.010);
1962
1963 TString text;
1964 text += (int)(range*.3);
1965 text += "mm";
1966
1967 TText newtxt2;
1968 newtxt2.SetTextSize(0.04);
1969 newtxt2.PaintText(-range*.95, -range*.95, text);
1970
1971 text = MString::Format("%.2f#circ", range*(0.99-0.60)*fGeomCam->GetConvMm2Deg());
1972 text = text.Strip(TString::kLeading);
1973
1974 TLatex latex;
1975 latex.PaintLatex(-range*.95, -range*.85, 0, 0.04, text);
1976 }
1977
1978 if (!TestBit(kNoLegend))
1979 {
1980 TPaveStats *stats = GetStatisticBox();
1981
1982 const Float_t hndc = 0.88 - (stats ? stats->GetY1NDC() : 1);
1983 const Float_t H = (0.80-hndc)*range;
1984 const Float_t offset = hndc*range;
1985
1986 const Int_t ncol = GetContour()==0 ? 50 : GetContour();
1987
1988 const Float_t h = 2./ncol;
1989 const Float_t w = range/sqrt((float)(fNcells-2));
1990
1991 TBox newbox;
1992 TText newtxt;
1993 newtxt.SetTextSize(0.03);
1994 newtxt.SetTextAlign(12);
1995#if ROOT_VERSION_CODE > ROOT_VERSION(3,01,06)
1996 newtxt.SetBit(/*kNoContextMenu|*/kCannotPick);
1997 newbox.SetBit(/*kNoContextMenu|*/kCannotPick);
1998#endif
1999
2000 const Float_t step = (islog && min>0 ? log10(max/min) : max-min);
2001 const Int_t firsts = step/48*3 < 1e-8 ? 8 : (Int_t)floor(log10(step/48*3));
2002 const TString opt = MString::Format("%%.%if", firsts>0 ? 0 : TMath::Abs(firsts));
2003/*
2004 for (Int_t i=0; i<ncol+1; i+=3)
2005 {
2006 Float_t val;
2007 if (islog && min>0)
2008 val = pow(10, step*i) * min;
2009 else
2010 val = min + step*i;
2011
2012 //const bool dispexp = max-min>1.5 && fabs(val)>0.1 && fabs(val)<1e6;
2013 newtxt.PaintText(range+1.5*w, H*(i*h-1)-offset, Form(opt, val));
2014 }
2015 */
2016 const MBinning bins(25, min, max, islog&&min>0?"log":"lin");
2017
2018 for (Int_t i=0; i<=25; i++)
2019 newtxt.PaintText(range+1.5*w, H*(i*ncol/25*h-1)-offset, MString::Format(opt, bins[i]));
2020
2021 for (Int_t i=0; i<ncol; i++)
2022 {
2023 newbox.SetFillColor(gStyle->GetColorPalette(i));
2024 newbox.PaintBox(range, H*(i*h-1)-offset, range+w, H*((i+1)*h-1)-offset);
2025 }
2026 }
2027}
2028
2029// ------------------------------------------------------------------------
2030//
2031// Save primitive as a C++ statement(s) on output stream out
2032//
2033void MHCamera::SavePrimitive(ostream &out, Option_t *opt)
2034{
2035 gLog << err << "MHCamera::SavePrimitive: Must be rewritten!" << endl;
2036 /*
2037 if (!gROOT->ClassSaved(TCanvas::Class()))
2038 fDrawingPad->SavePrimitive(out, opt);
2039
2040 out << " " << fDrawingPad->GetName() << "->SetWindowSize(";
2041 out << fDrawingPad->GetWw() << "," << fDrawingPad->GetWh() << ");" << endl;
2042 */
2043}
2044
2045void MHCamera::SavePrimitive(ofstream &out, Option_t *)
2046{
2047 MHCamera::SavePrimitive(static_cast<ostream&>(out), "");
2048}
2049
2050// ------------------------------------------------------------------------
2051//
2052// compute the distance of a point (px,py) to the Camera
2053// this functions needed for graphical primitives, that
2054// means without this function you are not able to interact
2055// with the graphical primitive with the mouse!!!
2056//
2057// All calcutations are done in pixel coordinates
2058//
2059Int_t MHCamera::DistancetoPrimitive(Int_t px, Int_t py)
2060{
2061 if (fNcells<=1)
2062 return 999999;
2063
2064 TPaveStats *box = (TPaveStats*)gPad->GetPrimitive("stats");
2065 if (box)
2066 {
2067 const Double_t w = box->GetY2NDC()-box->GetY1NDC();
2068 box->SetX1NDC(gStyle->GetStatX()-gStyle->GetStatW());
2069 box->SetY1NDC(gStyle->GetStatY()-w);
2070 box->SetX2NDC(gStyle->GetStatX());
2071 box->SetY2NDC(gStyle->GetStatY());
2072 }
2073
2074 if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
2075 return TH1D::DistancetoPrimitive(px, py);
2076
2077 const Bool_t issame = TString(GetDrawOption()).Contains("same", TString::kIgnoreCase);
2078
2079 const Float_t maxr = (1-fGeomCam->GetConvMm2Deg())*fGeomCam->GetMaxRadius()/2;
2080 const Float_t conv = !issame ||
2081 gPad->GetX1()<-maxr || gPad->GetY1()<-maxr ||
2082 gPad->GetX2()> maxr || gPad->GetY2()>maxr ? 1 : fGeomCam->GetConvMm2Deg();
2083
2084 if (GetPixelIndex(px, py, conv)>=0)
2085 return 0;
2086
2087 if (!box)
2088 return 999999;
2089
2090 const Int_t dist = box->DistancetoPrimitive(px, py);
2091 if (dist > TPad::GetMaxPickDistance())
2092 return 999999;
2093
2094 gPad->SetSelected(box);
2095 return dist;
2096}
2097
2098// ------------------------------------------------------------------------
2099//
2100//
2101Int_t MHCamera::GetPixelIndex(Int_t px, Int_t py, Float_t conv) const
2102{
2103 if (fNcells<=1)
2104 return -1;
2105
2106 for (Int_t i=0; i<fNcells-2; i++)
2107 if ((*fGeomCam)[i].DistancetoPrimitive(TMath::Nint(px*conv), TMath::Nint(py*conv))<=0)
2108 return i;
2109
2110 return -1;
2111}
2112
2113// ------------------------------------------------------------------------
2114//
2115// Returns string containing info about the object at position (px,py).
2116// Returned string will be re-used (lock in MT environment).
2117//
2118char *MHCamera::GetObjectInfo(Int_t px, Int_t py) const
2119{
2120 if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
2121#if ROOT_VERSION_CODE < ROOT_VERSION(5,22,00)
2122 return MH::GetObjectInfoH(px, py, *this);
2123#else
2124 return TH1D::GetObjectInfo(px, py);
2125#endif
2126
2127 static char info[128];
2128
2129 const Int_t idx=GetPixelIndex(px, py);
2130
2131 if (idx<0)
2132 return TObject::GetObjectInfo(px, py);
2133
2134 sprintf(info, "Software Pixel Idx: %d (Hardware Id=%d) c=%.1f <%s>",
2135 idx, idx+1, GetBinContent(idx+1), IsUsed(idx)?"on":"off");
2136 return info;
2137}
2138
2139// ------------------------------------------------------------------------
2140//
2141// Add a MCamEvent which should be displayed when the user clicks on a
2142// pixel.
2143// Warning: The object MUST inherit from TObject AND MCamEvent
2144//
2145void MHCamera::AddNotify(TObject *obj)
2146{
2147 // Make sure, that the object derives from MCamEvent!
2148 MCamEvent *evt = dynamic_cast<MCamEvent*>(obj);
2149 if (!evt)
2150 {
2151 gLog << err << "ERROR: MHCamera::AddNotify - TObject doesn't inherit from MCamEvent... ignored." << endl;
2152 return;
2153 }
2154
2155 // Make sure, that it is deleted from the list too, if the obj is deleted
2156 obj->SetBit(kMustCleanup);
2157
2158 // Add object to list
2159 fNotify->Add(obj);
2160}
2161
2162// ------------------------------------------------------------------------
2163//
2164// Execute a mouse event on the camera
2165//
2166void MHCamera::ExecuteEvent(Int_t event, Int_t px, Int_t py)
2167{
2168 if (TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase))
2169 {
2170 TH1D::ExecuteEvent(event, px, py);
2171 return;
2172 }
2173 //if (event==kMouseMotion && fStatusBar)
2174 // fStatusBar->SetText(GetObjectInfo(px, py), 0);
2175 if (event!=kButton1Down)
2176 return;
2177
2178 const Int_t idx = GetPixelIndex(px, py);
2179 if (idx<0)
2180 return;
2181
2182 gLog << all << GetTitle() << " <" << GetName() << ">" << dec << endl;
2183 gLog << "Software Pixel Idx: " << idx << endl;
2184 gLog << "Hardware Pixel Id: " << idx+1 << endl;
2185 gLog << "Contents: " << GetBinContent(idx+1);
2186 if (GetBinError(idx+1)>0)
2187 gLog << " +/- " << GetBinError(idx+1);
2188 gLog << " <" << (IsUsed(idx)?"on":"off") << "> n=" << fBinEntries[idx+1] << endl;
2189
2190 if (fNotify && fNotify->GetSize()>0)
2191 {
2192 // FIXME: Is there a simpler and more convinient way?
2193
2194 // The name which is created here depends on the instance of
2195 // MHCamera and on the pad on which it is drawn --> The name
2196 // is unique. For ExecuteEvent gPad is always correctly set.
2197 const TString name = MString::Format("%p;%p;PixelContent", this, gPad);
2198
2199 TCanvas *old = (TCanvas*)gROOT->GetListOfCanvases()->FindObject(name);
2200 if (old)
2201 old->cd();
2202 else
2203 new TCanvas(name);
2204
2205 /*
2206 TIter Next(gPad->GetListOfPrimitives());
2207 TObject *o;
2208 while (o=Next()) cout << o << ": " << o->GetName() << " " << o->IsA()->GetName() << endl;
2209 */
2210
2211 // FIXME: Make sure, that the old histograms are really deleted.
2212 // Are they already deleted?
2213
2214 // The dynamic_cast is necessary here: We cannot use ForEach
2215 TIter Next(fNotify);
2216 MCamEvent *evt;
2217 while ((evt=dynamic_cast<MCamEvent*>(Next())))
2218 evt->DrawPixelContent(idx);
2219
2220 gPad->Modified();
2221 gPad->Update();
2222 }
2223}
2224
2225UInt_t MHCamera::GetNumPixels() const
2226{
2227 return fGeomCam ? fGeomCam->GetNumPixels() : 0;
2228}
2229
2230// --------------------------------------------------------------------------
2231//
2232// Draw a projection of MHCamera onto the y-axis values. Depending on the
2233// variable fit, the following fits are performed:
2234//
2235// 0: No fit, simply draw the projection
2236// 1: Single Gauss (for distributions flat-fielded over the whole camera)
2237// 2: Double Gauss (for distributions different for inner and outer pixels)
2238// 3: Triple Gauss (for distributions with inner, outer pixels and outliers)
2239// 4: flat (for the probability distributions)
2240// (1-4:) Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are
2241// drawn separately, for inner and outer pixels.
2242// 5: Fit Inner and Outer pixels separately by a single Gaussian
2243// (only for MAGIC cameras)
2244// 6: Fit Inner and Outer pixels separately by a single Gaussian and display
2245// additionally the two camera halfs separately (for MAGIC camera)
2246// 7: Single Gauss with TLegend to show the meaning of the colours
2247//
2248void MHCamera::DrawProjection(Int_t fit) const
2249{
2250 if (fit==5 || fit==6)
2251 {
2252 const UInt_t n = fGeomCam->GetNumAreas();
2253
2254 TVirtualPad *pad = gPad;
2255 pad->Divide(n, 1, 1e-5, 1e-5);;
2256
2257 for (UInt_t i=0; i<n; i++)
2258 {
2259 pad->cd(i+1);
2260 gPad->SetBorderMode(0);
2261 gPad->SetRightMargin(0.025);
2262 gPad->SetTicks();
2263
2264 TH1D &h = *ProjectionS(TArrayI(), TArrayI(1, (Int_t*)&i),
2265 MString::Format("%s_%d", GetName(), i));
2266 h.SetTitle(MString::Format("%s %d", GetTitle(), i));
2267 h.SetDirectory(NULL);
2268 h.SetBit(kCanDelete);
2269 h.Draw();
2270
2271 TAxis *xaxe = h.GetXaxis();
2272 TAxis *yaxe = h.GetYaxis();
2273
2274 xaxe->CenterTitle();
2275 yaxe->CenterTitle();
2276 xaxe->SetTitleSize(0.06);
2277 yaxe->SetTitleSize(0.06);
2278 xaxe->SetTitleOffset(0.8);
2279 yaxe->SetTitleOffset(0.85);
2280 xaxe->SetLabelSize(0.05);
2281 yaxe->SetLabelSize(0.05);
2282 if (i>0)
2283 yaxe->SetTitle("");
2284
2285 if (fit==5)
2286 continue;
2287
2288 h.Fit("gaus", "Q");
2289
2290 TF1 *f = h.GetFunction("gaus");
2291 if (f)
2292 {
2293 f->SetLineWidth(2);
2294 f->SetLineColor(kBlue);
2295 }
2296 }
2297 return;
2298 }
2299
2300 TH1D *obj2 = (TH1D*)Projection(GetName());
2301 obj2->SetDirectory(0);
2302 obj2->Draw();
2303 obj2->SetBit(kCanDelete);
2304
2305 if (fit==0)
2306 return;
2307
2308 const Double_t xmin = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
2309 const Double_t xmax = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
2310 const Double_t integ = obj2->Integral("width")/2.5;
2311 const Double_t max = obj2->GetMaximum();
2312 const Double_t mean = obj2->GetMean();
2313 const Double_t rms = obj2->GetRMS();
2314 const Double_t width = xmax-xmin;
2315
2316 const char *dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
2317 "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
2318
2319 const char *tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
2320 "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"
2321 "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
2322
2323 TF1 *f=0;
2324 switch (fit)
2325 {
2326 case 1: // Single gauss
2327 f = new TF1("sgaus", "gaus", xmin, xmax);
2328 f->SetLineColor(kBlue);
2329 f->SetBit(kCanDelete);
2330 f->SetParNames("Max", "#mu", "#sigma");
2331 f->SetParameters(max, mean, rms);
2332 f->SetParLimits(0, 0, max*2);
2333 f->SetParLimits(1, xmin, xmax);
2334 f->SetParLimits(2, 0, width/1.5);
2335
2336 obj2->Fit(f, "QLR");
2337 break;
2338
2339 case 2: // Double gauss
2340 f = new TF1("dgaus", dgausformula, xmin, xmax);
2341 f->SetLineColor(kBlue);
2342 f->SetBit(kCanDelete);
2343 f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
2344 f->SetParameters(integ, (xmin+mean)/2., width/4.,
2345 integ/width/2., (xmax+mean)/2., width/4.);
2346 // The left-sided Gauss
2347 f->SetParLimits(0, integ-1.5, integ+1.5);
2348 f->SetParLimits(1, xmin+(width/10.), mean);
2349 f->SetParLimits(2, 0, width/2.);
2350 // The right-sided Gauss
2351 f->SetParLimits(3, 0, integ);
2352 f->SetParLimits(4, mean, xmax-(width/10.));
2353 f->SetParLimits(5, 0, width/2.);
2354 obj2->Fit(f,"QLRM");
2355 break;
2356
2357 case 3: // Triple gauss
2358 f = new TF1("tgaus", tgausformula, xmin,xmax);
2359 f->SetLineColor(kBlue);
2360 f->SetBit(kCanDelete);
2361 f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
2362 "A_{2}","#mu_{2}","#sigma_{2}",
2363 "A_{3}","#mu_{3}","#sigma_{3}");
2364 f->SetParameters(integ, (xmin+mean)/2, width/4.,
2365 integ/width/3., (xmax+mean)/2., width/4.,
2366 integ/width/3., mean,width/2.);
2367 // The left-sided Gauss
2368 f->SetParLimits(0, integ-1.5, integ+1.5);
2369 f->SetParLimits(1, xmin+(width/10.), mean);
2370 f->SetParLimits(2, width/15., width/2.);
2371 // The right-sided Gauss
2372 f->SetParLimits(3, 0., integ);
2373 f->SetParLimits(4, mean, xmax-(width/10.));
2374 f->SetParLimits(5, width/15., width/2.);
2375 // The Gauss describing the outliers
2376 f->SetParLimits(6, 0., integ);
2377 f->SetParLimits(7, xmin, xmax);
2378 f->SetParLimits(8, width/4., width/1.5);
2379 obj2->Fit(f,"QLRM");
2380 break;
2381
2382 case 4:
2383 obj2->Fit("pol0", "Q");
2384 if (obj2->GetFunction("pol0"))
2385 obj2->GetFunction("pol0")->SetLineColor(kBlue);
2386 break;
2387
2388 case 9:
2389 break;
2390
2391 default:
2392 obj2->Fit("gaus", "Q");
2393 if (obj2->GetFunction("gaus"))
2394 obj2->GetFunction("gaus")->SetLineColor(kBlue);
2395 break;
2396 }
2397}
2398
2399// --------------------------------------------------------------------------
2400//
2401// Draw a projection of MHCamera vs. the radius from the central pixel.
2402//
2403// The inner and outer pixels are drawn separately, both fitted by a polynomial
2404// of grade 1.
2405//
2406void MHCamera::DrawRadialProfile() const
2407{
2408 TProfile *obj2 = (TProfile*)RadialProfile(GetName());
2409 obj2->SetDirectory(0);
2410 obj2->Draw();
2411 obj2->SetBit(kCanDelete);
2412/*
2413 if (GetGeomCam().InheritsFrom("MGeomCamMagic"))
2414 {
2415 TArrayI s0(6);
2416 s0[0] = 1;
2417 s0[1] = 2;
2418 s0[2] = 3;
2419 s0[3] = 4;
2420 s0[4] = 5;
2421 s0[5] = 6;
2422
2423 TArrayI inner(1);
2424 inner[0] = 0;
2425
2426 TArrayI outer(1);
2427 outer[0] = 1;
2428
2429 // Just to get the right (maximum) binning
2430 TProfile *half[2];
2431 half[0] = RadialProfileS(s0, inner, MString::Format("%sInner",GetName()));
2432 half[1] = RadialProfileS(s0, outer, MString::Format("%sOuter",GetName()));
2433
2434 for (Int_t i=0; i<2; i++)
2435 {
2436 Double_t min = GetGeomCam().GetMinRadius(i);
2437 Double_t max = GetGeomCam().GetMaxRadius(i);
2438
2439 half[i]->SetLineColor(kRed+i);
2440 half[i]->SetDirectory(0);
2441 half[i]->SetBit(kCanDelete);
2442 half[i]->Draw("same");
2443 half[i]->Fit("pol1","Q","",min,max);
2444 half[i]->GetFunction("pol1")->SetLineColor(kRed+i);
2445 half[i]->GetFunction("pol1")->SetLineWidth(1);
2446 }
2447 }
2448*/
2449}
2450
2451// --------------------------------------------------------------------------
2452//
2453// Draw a projection of MHCamera vs. the azimuth angle inside the camera.
2454//
2455// The inner and outer pixels are drawn separately.
2456// The general azimuth profile is fitted by a straight line
2457//
2458void MHCamera::DrawAzimuthProfile() const
2459{
2460 TProfile *obj2 = (TProfile*)AzimuthProfile(GetName());
2461 obj2->SetDirectory(0);
2462 obj2->Draw();
2463 obj2->SetBit(kCanDelete);
2464 obj2->Fit("pol0","Q","");
2465 if (obj2->GetFunction("pol0"))
2466 obj2->GetFunction("pol0")->SetLineWidth(1);
2467/*
2468 if (GetGeomCam().InheritsFrom("MGeomCamMagic"))
2469 {
2470 TArrayI inner(1);
2471 inner[0] = 0;
2472
2473 TArrayI outer(1);
2474 outer[0] = 1;
2475
2476 // Just to get the right (maximum) binning
2477 TProfile *half[2];
2478 half[0] = AzimuthProfileA(inner, MString::Format("%sInner",GetName()));
2479 half[1] = AzimuthProfileA(outer, MString::Format("%sOuter",GetName()));
2480
2481 for (Int_t i=0; i<2; i++)
2482 {
2483 half[i]->SetLineColor(kRed+i);
2484 half[i]->SetDirectory(0);
2485 half[i]->SetBit(kCanDelete);
2486 half[i]->SetMarkerSize(0.5);
2487 half[i]->Draw("same");
2488 }
2489 }
2490*/
2491}
2492
2493// --------------------------------------------------------------------------
2494//
2495// Draw the MHCamera into the MStatusDisplay:
2496//
2497// 1) Draw it as histogram (MHCamera::DrawCopy("hist")
2498// 2) Draw it as a camera, with MHCamera::SetPrettyPalette() set.
2499// 3) If "rad" is not zero, draw its values vs. the radius from the camera center.
2500// (DrawRadialProfile())
2501// 4) Depending on the variable "fit", draw the values projection on the y-axis
2502// (DrawProjection()):
2503// 0: don't draw
2504// 1: Draw fit to Single Gauss (for distributions flat-fielded over the whole camera)
2505// 2: Draw and fit to Double Gauss (for distributions different for inner and outer pixels)
2506// 3: Draw and fit to Triple Gauss (for distributions with inner, outer pixels and outliers)
2507// 4: Draw and fit to Polynomial grade 0: (for the probability distributions)
2508// >4: Draw and don;t fit.
2509//
2510void MHCamera::CamDraw(TCanvas &c, const Int_t x, const Int_t y,
2511 const Int_t fit, const Int_t rad, const Int_t azi,
2512 TObject *notify)
2513{
2514 c.cd(x);
2515 gPad->SetBorderMode(0);
2516 gPad->SetRightMargin(0.02);
2517 gPad->SetTicks();
2518 MHCamera *obj1=(MHCamera*)DrawCopy("hist");
2519 obj1->SetDirectory(NULL);
2520 obj1->SetStats(kFALSE);
2521
2522 if (notify)
2523 obj1->AddNotify(notify);
2524
2525 c.cd(x+y);
2526 gPad->SetBorderMode(0);
2527 obj1->SetPrettyPalette();
2528 obj1->Draw();
2529
2530 Int_t cnt = 2;
2531
2532 if (rad)
2533 {
2534 c.cd(x+2*y);
2535 gPad->SetBorderMode(0);
2536 gPad->SetTicks();
2537 DrawRadialProfile();
2538 cnt++;
2539 }
2540
2541 if (azi)
2542 {
2543 c.cd(x+cnt*y);
2544 gPad->SetBorderMode(0);
2545 gPad->SetTicks();
2546 DrawAzimuthProfile();
2547 cnt++;
2548 }
2549
2550 if (fit<0)
2551 return;
2552
2553 c.cd(x + cnt*y);
2554 gPad->SetBorderMode(0);
2555 gPad->SetRightMargin(0.025);
2556 gPad->SetTicks();
2557 DrawProjection(fit);
2558}
Note: See TracBrowser for help on using the repository browser.