source: trunk/MagicSoft/Mars/mhist/MHCamera.cc@ 9369

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