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

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