source: releases/Mars.2014.05.26/mhist/MHCamera.cc@ 19403

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