source: trunk/Mars/mhflux/MHAlpha.cc@ 19187

Last change on this file since 19187 was 18957, checked in by tbretz, 7 years ago
Improved compiler warnings.
File size: 32.2 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, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MHAlpha
28//
29// Create a single Alpha-Plot. The alpha-plot is fitted online. You can
30// check the result when it is filles in the MStatusDisplay
31// For more information see MHFalseSource::FitSignificance
32//
33// For convinience (fit) the output significance is stored in a
34// container in the parlisrt
35//
36// PRELIMINARY!
37//
38//
39// ToDo:
40// =====
41//
42// - Replace time-binning by histogram (which is enhanced bin-wise)
43//
44//
45//////////////////////////////////////////////////////////////////////////////
46#include "MHAlpha.h"
47
48#include <TStyle.h>
49#include <TLatex.h>
50#include <TCanvas.h>
51#include <TPaveStats.h>
52
53#include "MHillas.h"
54#include "MHillasSrc.h"
55#include "MTime.h"
56#include "MObservatory.h"
57#include "MPointingPos.h"
58#include "MAstroSky2Local.h"
59#include "MStatusDisplay.h"
60#include "MParameters.h"
61#include "MHMatrix.h"
62
63#include "MString.h"
64#include "MBinning.h"
65#include "MParList.h"
66
67#include "MLog.h"
68#include "MLogManip.h"
69
70ClassImp(MHAlpha);
71
72using namespace std;
73
74// --------------------------------------------------------------------------
75//
76// Default Constructor
77//
78MHAlpha::MHAlpha(const char *name, const char *title)
79 : fNameParameter("MHillasSrc"), fParameter(0),
80 fOffData(0), fResult(0), fSigma(0), fEnergy(0), fBin(0),
81 fPointPos(0), fTimeEffOn(0), fTime(0), fNumTimeBins(10),
82 fHillas(0), fMatrix(0), fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE),
83 fSkipHistEnergy(kFALSE), fForceUsingSize(kFALSE)
84{
85 //
86 // set the name and title of this object
87 //
88 fName = name ? name : "MHAlpha";
89 fTitle = title ? title : "Alpha plot";
90
91 fHist.SetName("Alpha");
92 fHist.SetTitle("Alpha");
93 fHist.SetXTitle("\\Theta [deg]");
94 //fHist.SetYTitle("E_{est} [GeV]");
95 fHist.SetZTitle("|\\alpha| [\\circ]");
96 fHist.SetDirectory(NULL);
97 fHist.UseCurrentStyle();
98
99 // Main histogram
100 fHistTime.SetName("Alpha");
101 fHistTime.SetXTitle("|\\alpha| [\\circ]");
102 fHistTime.SetYTitle("Counts");
103 fHistTime.UseCurrentStyle();
104 fHistTime.SetDirectory(NULL);
105
106 fHEnergy.SetName("Excess");
107 //fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
108 //fHEnergy.SetXTitle("E_{est} [GeV]");
109 fHEnergy.SetYTitle("N_{exc}");
110 fHEnergy.SetDirectory(NULL);
111 fHEnergy.UseCurrentStyle();
112
113 fHTheta.SetName("ExcessTheta");
114 fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
115 fHTheta.SetXTitle("\\Theta [\\circ]");
116 fHTheta.SetYTitle("N_{exc}");
117 fHTheta.SetDirectory(NULL);
118 fHTheta.UseCurrentStyle();
119 fHTheta.SetMinimum(0);
120
121 // effective on time versus time
122 fHTime.SetName("ExcessTime");
123 fHTime.SetTitle(" N_{exc} vs. Time ");
124 fHTime.SetXTitle("Time");
125 fHTime.SetYTitle("N_{exc} [s]");
126 fHTime.UseCurrentStyle();
127 fHTime.SetDirectory(NULL);
128 fHTime.GetYaxis()->SetTitleOffset(1.2);
129 fHTime.GetXaxis()->SetLabelSize(0.033);
130 fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
131 fHTime.GetXaxis()->SetTimeDisplay(1);
132 fHTime.SetMinimum(0);
133 fHTime.Sumw2();
134
135 MBinning binsa, binse, binst;
136 binsa.SetEdges(18, 0, 90);
137 binse.SetEdgesLog(15, 10, 100000);
138 binst.SetEdgesASin(67, -0.005, 0.665);
139 binse.Apply(fHEnergy);
140 binst.Apply(fHTheta);
141 binsa.Apply(fHistTime);
142
143 MH::SetBinning(fHist, binst, binse, binsa);
144}
145
146Float_t MHAlpha::FitEnergyBins(Bool_t paint)
147{
148 // Do not store the 'final' result in fFit
149 MAlphaFitter fit(fFit);
150
151 const Int_t n = fHist.GetNbinsY();
152
153 fHEnergy.SetEntries(0);
154
155 Float_t mean = 0;
156 Int_t num = 0;
157 for (int i=1; i<=n; i++)
158 {
159 if (!fit.FitEnergy(fHist, fOffData, i))
160 continue;
161
162 // FIXME: Calculate UL!
163 if (fit.GetEventsExcess()<=0)
164 continue;
165
166 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
167 fHEnergy.SetBinError(i, fit.GetErrorExcess());
168
169 mean += fit.GetEventsExcess()*fit.GetEventsExcess()/fit.GetErrorExcess()/fit.GetErrorExcess();
170 num++;
171 }
172 return TMath::Sqrt(mean)/num;
173}
174
175void MHAlpha::FitThetaBins(Bool_t paint)
176{
177 // Do not store the 'final' result in fFit
178 MAlphaFitter fit(fFit);
179
180 const Int_t n = fHist.GetNbinsX();
181
182 fHTheta.SetEntries(0);
183
184 for (int i=1; i<=n; i++)
185 {
186 if (!fit.FitTheta(fHist, fOffData, i))
187 continue;
188
189 // FIXME: Calculate UL!
190 if (fit.GetEventsExcess()<=0)
191 continue;
192
193 fHTheta.SetBinContent(i, fit.GetEventsExcess());
194 fHTheta.SetBinError(i, fit.GetErrorExcess());
195 }
196}
197
198// --------------------------------------------------------------------------
199//
200// Return the value from fParemeter which should be filled into the plots
201//
202Double_t MHAlpha::GetVal() const
203{
204 return static_cast<const MHillasSrc*>(fParameter)->GetAlpha();
205}
206
207
208// --------------------------------------------------------------------------
209//
210// Store the pointer to the parameter container storing the plotted value
211// (here MHillasSrc) in fParameter.
212//
213// return whether it was found or not.
214//
215Bool_t MHAlpha::GetParameter(const MParList &pl)
216{
217 fParameter = (MParContainer*)pl.FindObject(fNameParameter, "MHillasSrc");
218 if (fParameter)
219 return kTRUE;
220
221 *fLog << err << fNameParameter << " [MHillasSrc] not found... abort." << endl;
222 return kFALSE;
223}
224
225Bool_t MHAlpha::SetupFill(const MParList *pl)
226{
227 fHist.Reset();
228 fHEnergy.Reset();
229 fHTheta.Reset();
230 fHTime.Reset();
231
232 const TString off(MString::Format("%sOff", fName.Data()));
233 if (fName!=off && fOffData==NULL)
234 {
235 const TString desc(MString::Format("%s [%s] found... using ", off.Data(), ClassName()));
236 MHAlpha *hoff = (MHAlpha*)pl->FindObject(off, ClassName());
237 if (!hoff)
238 *fLog << inf << "No " << desc << "current data only!" << endl;
239 else
240 {
241 *fLog << inf << desc << "on-off mode!" << endl;
242 SetOffData(*hoff);
243 }
244 }
245
246 if (!GetParameter(*pl))
247 return kFALSE;
248
249 fHillas = 0;
250 fEnergy = fForceUsingSize ? 0 : (MParameterD*)pl->FindObject("MEnergyEst", "MParameterD");
251 if (!fEnergy)
252 {
253 *fLog << warn << "MEnergyEst [MParameterD] not found... " << flush;
254
255 if (!fHillas)
256 fHillas = (MHillas*)pl->FindObject("MHillas");
257 if (fHillas)
258 *fLog << "using SIZE instead!" << endl;
259 else
260 *fLog << "ignored." << endl;
261
262 fHEnergy.SetTitle(" N_{exc} vs. Size ");
263 fHEnergy.SetXTitle("Size [phe]");
264 fHist.SetYTitle("Size [phe]");
265 }
266 else
267 {
268 fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
269 fHEnergy.SetXTitle("E_{est} [GeV]");
270 fHist.SetYTitle("E_{est} [GeV]");
271 }
272
273 fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
274 if (!fPointPos)
275 *fLog << warn << "MPointingPos not found... ignored." << endl;
276
277 fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
278 if (!fTimeEffOn)
279 *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
280 else
281 *fTimeEffOn = MTime(); // FIXME: How to do it different?
282
283 fTime = (MTime*)pl->FindObject("MTime");
284 if (!fTime)
285 *fLog << warn << "MTime not found... ignored." << endl;
286
287 fResult = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MinimizationValue");
288 if (!fResult)
289 return kFALSE;
290 fSigma = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "GaussSigma");
291 if (!fSigma)
292 return kFALSE;
293 fBin = (MParameterI*)const_cast<MParList*>(pl)->FindCreateObj("MParameterI", "Bin");
294 if (!fBin)
295 return kFALSE;
296
297 //fExcess = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MExcess");
298 //if (!fExcess)
299 // return kFALSE;
300
301 fLastTime = MTime();
302 fNumRebin = fNumTimeBins;
303
304 MBinning binst, binse, binsa;
305 binst.SetEdges(fOffData ? *fOffData : fHist, 'x');
306 binse.SetEdges(fOffData ? *fOffData : fHist, 'y');
307 binsa.SetEdges(fOffData ? *fOffData : fHist, 'z');
308 if (!fOffData)
309 {
310 if (!fPointPos)
311 binst.SetEdges(1, 0, 60);
312 else
313 binst.SetEdges(*pl, "BinningTheta");
314
315 if (!fEnergy && !fHillas)
316 binse.SetEdges(1, 10, 100000);
317 else
318 if (fEnergy)
319 binse.SetEdges(*pl, "BinningEnergyEst");
320 else
321 binse.SetEdges(*pl, "BinningSize");
322
323 binsa.SetEdges(*pl, MString::Format("Binning%s", ClassName()+2));
324 }
325 else
326 {
327 fHEnergy.SetTitle(fOffData->GetTitle());
328 fHEnergy.SetXTitle(fOffData->GetYaxis()->GetTitle());
329 fHist.SetYTitle(fOffData->GetYaxis()->GetTitle());
330 }
331
332 binse.Apply(fHEnergy);
333 binst.Apply(fHTheta);
334 binsa.Apply(fHistTime);
335 MH::SetBinning(fHist, binst, binse, binsa);
336
337 MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
338 if (!fit)
339 *fLog << warn << "MAlphaFitter not found... ignored." << endl;
340 else
341 fit->Copy(fFit);
342
343 *fLog << inf;
344 fFit.Print();
345
346 return kTRUE;
347}
348
349void MHAlpha::InitAlphaTime(const MTime &t)
350{
351 //
352 // If this is the first call we have to initialize the time-histogram
353 //
354 const MBinning bins(1, t.GetAxisTime()-60, t.GetAxisTime());
355 bins.Apply(fHTime);
356
357 fLastTime=t;
358}
359
360void MHAlpha::UpdateAlphaTime(Bool_t final)
361{
362 if (!fTimeEffOn)
363 return;
364
365 if (!final)
366 {
367 if (fLastTime==MTime() && *fTimeEffOn!=MTime())
368 {
369 InitAlphaTime(*fTimeEffOn);
370 return;
371 }
372
373 if (fLastTime!=*fTimeEffOn)
374 {
375 fLastTime=*fTimeEffOn;
376 fNumRebin--;
377 }
378
379 if (fNumRebin>0)
380 return;
381 }
382
383 // Check new 'last time'
384 MTime *time = final ? fTime : fTimeEffOn;
385
386 if (time->GetAxisTime()<=fHTime.GetXaxis()->GetXmax())
387 {
388 *fLog << warn << "WARNING - New time-stamp " << *time << " lower" << endl;
389 *fLog << "than upper edge of histogram... skipped." << endl;
390 *fLog << "This should not happen. Maybe you started you eventloop with" << endl;
391 *fLog << "an already initialized time stamp MTimeEffectiveOnTime?" << endl;
392 fNumRebin++;
393 return;
394 }
395
396 // Fit time histogram
397 MAlphaFitter fit(fFit);
398
399 TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E") : 0;
400 const Bool_t rc = fit.ScaleAndFit(fHistTime, h);
401
402 if (h)
403 delete h;
404
405 if (!rc)
406 return;
407
408 // Reset Histogram
409 fHistTime.Reset();
410 fHistTime.SetEntries(0);
411
412 //
413 // Prepare Histogram
414 //
415 if (final)
416 time->Plus1ns();
417
418 // Enhance binning
419 MBinning bins;
420 bins.SetEdges(fHTime, 'x');
421 bins.AddEdge(time->GetAxisTime());
422 bins.Apply(fHTime);
423
424 //
425 // Fill histogram
426 //
427 // Get number of bins
428 const Int_t n = fHTime.GetNbinsX();
429
430 fHTime.SetBinContent(n, fit.GetEventsExcess());
431 fHTime.SetBinError(n, fit.GetErrorExcess());
432
433 //*fLog << inf3 << *fTimeEffOn << " (" << n << "): " << fit.GetEventsExcess() << endl;
434
435 fNumRebin = fNumTimeBins;
436}
437
438void MHAlpha::SetBin(Int_t ibin)
439{
440 // Is this necessary?
441 // Could be speed up up searching for it only once.
442 const Float_t max = fFit.GetSignalIntegralMax();
443 const Int_t bin0 = fHist.GetZaxis()->FindFixBin(max);
444
445 const Int_t nbinsx = fHist.GetNbinsX();
446 const Int_t nbinsy = fHist.GetNbinsY();
447 const Int_t nxy = (nbinsx+2)*(nbinsy+2);
448
449 const Int_t binz = ibin/nxy;
450
451 const Bool_t issignal = binz>0 && binz<bin0;
452
453 fBin->SetVal(issignal ? binz : -binz);
454}
455
456// --------------------------------------------------------------------------
457//
458// Fill the histogram. For details see the code or the class description
459//
460Int_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
461{
462 Double_t alpha, energy, theta;
463 Double_t size=-1;
464
465 if (fMatrix)
466 {
467 alpha = fMap[0]<0 ? GetVal() : (*fMatrix)[fMap[0]];
468 energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]];
469 size = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]];
470 //<0 ? 1000 : (*fMatrix)[fMap[1]];
471 theta = 0;
472
473 if (energy<0)
474 energy=size;
475 if (size<0)
476 size=energy;
477
478 if (energy<0 && size<0)
479 energy = size = 1000;
480 }
481 else
482 {
483 alpha = GetVal();
484
485 if (fHillas)
486 size = fHillas->GetSize();
487 energy = fEnergy ? fEnergy->GetVal() : (fHillas?fHillas->GetSize():1000);
488 theta = fPointPos ? fPointPos->GetZd() : 0;
489 }
490
491 // enhance histogram if necessary
492 UpdateAlphaTime();
493
494 // Fill histograms
495 const Int_t ibin = fHist.Fill(theta, energy, TMath::Abs(alpha), w);
496 SetBin(ibin);
497
498 if (!fSkipHistTime)
499 fHistTime.Fill(TMath::Abs(alpha), w);
500
501 return kTRUE;
502}
503
504// --------------------------------------------------------------------------
505//
506// Paint the integral and the error on top of the histogram
507//
508void MHAlpha::PaintText(Double_t val, Double_t error) const
509{
510 TLatex text(0.45, 0.94, MString::Format("N_{exc} = %.1f \\pm %.1f", val, error));
511 text.SetBit(TLatex::kTextNDC);
512 text.SetTextSize(0.04);
513 text.Paint();
514}
515
516// --------------------------------------------------------------------------
517//
518// Paint the integral and the error on top of the histogram
519//
520void MHAlpha::PaintText(const TH1D &h) const
521{
522 Double_t sumv = 0;
523 Double_t sume = 0;
524
525 for (int i=0; i<h.GetNbinsX(); i++)
526 {
527 sumv += h.GetBinContent(i+1);
528 sume += h.GetBinError(i+1);
529 }
530 PaintText(sumv, sume);
531}
532// --------------------------------------------------------------------------
533//
534// Update the projections
535//
536void MHAlpha::Paint(Option_t *opt)
537{
538 // Note: Any object cannot be added to a pad twice!
539 // The object is searched by gROOT->FindObject only in
540 // gPad itself!
541 TVirtualPad *padsave = gPad;
542
543 TH1D *h0=0;
544
545 TString o(opt);
546
547 if (o==(TString)"proj")
548 {
549 TPaveStats *st=0;
550 for (int x=0; x<4; x++)
551 {
552 TVirtualPad *p=padsave->GetPad(x+1);
553 if (!p || !(st = (TPaveStats*)p->GetPrimitive("stats")))
554 continue;
555
556 if (st->GetOptStat()==11)
557 continue;
558
559 const Double_t y1 = st->GetY1NDC();
560 const Double_t y2 = st->GetY2NDC();
561 const Double_t x1 = st->GetX1NDC();
562 const Double_t x2 = st->GetX2NDC();
563
564 st->SetY1NDC((y2-y1)/3+y1);
565 st->SetX1NDC((x2-x1)/3+x1);
566 st->SetOptStat(11);
567 }
568
569 padsave->cd(1);
570
571 TH1D *hon = (TH1D*)gPad->FindObject("Proj");
572 if (hon)
573 {
574 TH1D *dum = fHist.ProjectionZ("dumab", 0, fHist.GetNbinsX()+1, 0, fHist.GetNbinsY()+1);
575 dum->SetDirectory(0);
576 hon->Reset();
577 hon->Add(dum);
578 delete dum;
579
580 if (fOffData)
581 {
582 TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
583 if (hoff)
584 {
585 TH1D *dum2 = fOffData->ProjectionZ("dumxy", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1);
586 dum2->SetDirectory(0);
587 hoff->Reset();
588 hoff->Add(dum2);
589 delete dum2;
590
591 const Double_t alpha = fFit.Scale(*hoff, *hon);
592
593 hon->SetMaximum();
594 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
595
596 // BE CARFEULL: This is a really weird workaround!
597 hoff->SetMaximum(alpha);
598
599 // For some reason the line-color is resetted
600 hoff->SetLineColor(kRed);
601
602 if ((h0=(TH1D*)gPad->FindObject("ProjOnOff")))
603 {
604 h0->Reset();
605 h0->Add(hoff, hon, -1);
606 const Float_t min = h0->GetMinimum()*1.05;
607 hon->SetMinimum(min<0 ? min : 0);
608 }
609
610 }
611 }
612 else
613 hon->SetMinimum(0);
614 }
615 FitEnergyBins();
616 FitThetaBins();
617 }
618
619 if (o==(TString)"variable")
620 if ((h0 = (TH1D*)gPad->FindObject("Proj")))
621 {
622 // Check whether we are dealing with on-off analysis
623 TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
624
625 // BE CARFEULL: This is a really weird workaround!
626 const Double_t scale = !hoff || hoff->GetMaximum()<0 ? 1 : hoff->GetMaximum();
627
628 // Do not store the 'final' result in fFit
629 MAlphaFitter fit(fFit);
630 fit.Fit(*h0, hoff, scale, kTRUE);
631 fit.PaintResult();
632 }
633
634 if (o==(TString)"time")
635 PaintText(fHTime);
636
637 if (o==(TString)"theta")
638 {
639 TH1 *h = (TH1*)gPad->FindObject(MString::Format("%s_x", fHist.GetName()));
640 if (h)
641 {
642 TH1D *h2 = (TH1D*)fHist.Project3D("dum_x");
643 h2->SetDirectory(0);
644 h2->Scale(fHTheta.Integral()/h2->Integral());
645 h->Reset();
646 h->Add(h2);
647 delete h2;
648 }
649 PaintText(fHTheta);
650 }
651
652 if (o==(TString)"energy")
653 {
654 TH1 *h = (TH1*)gPad->FindObject(MString::Format("%s_y", fHist.GetName()));
655 if (h)
656 {
657 TH1D *h2 = (TH1D*)fHist.Project3D("dum_y");
658 h2->SetDirectory(0);
659 h2->Scale(fHEnergy.Integral()/h2->Integral());
660 h->Reset();
661 h->Add(h2);
662 delete h2;
663 }
664 PaintText(fHEnergy);
665
666 if (fHEnergy.GetMaximum()>1)
667 {
668 gPad->SetLogx();
669 gPad->SetLogy();
670 }
671 }
672
673 gPad = padsave;
674}
675
676// --------------------------------------------------------------------------
677//
678// Draw the histogram
679//
680void MHAlpha::Draw(Option_t *opt)
681{
682 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
683
684 /*
685 if (TString(opt).Contains("sizebins", TString::kIgnoreCase))
686 {
687 AppendPad("sizebins");
688 return;
689 }
690 */
691
692 // Do the projection before painting the histograms into
693 // the individual pads
694 AppendPad("proj");
695
696 pad->SetBorderMode(0);
697 pad->Divide(2,2);
698
699 TH1D *h=0;
700
701 pad->cd(1);
702 gPad->SetBorderMode(0);
703
704 h = fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, 0, fHist.GetNbinsY()+1, "E");
705 h->SetBit(TH1::kNoTitle);
706 h->SetStats(kTRUE);
707 h->SetXTitle(fHist.GetZaxis()->GetTitle());
708 h->SetYTitle("Counts");
709 h->SetDirectory(NULL);
710 h->SetMarkerStyle(0);
711 h->SetBit(kCanDelete);
712 h->Draw("");
713
714 if (fOffData)
715 {
716 // To get green on-data
717 //h->SetMarkerColor(kGreen);
718 //h->SetLineColor(kGreen);
719
720 h = fOffData->ProjectionZ("ProjOff", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E");
721 h->SetBit(TH1::kNoTitle);
722 h->SetXTitle(fHist.GetZaxis()->GetTitle());
723 h->SetYTitle("Counts");
724 h->SetDirectory(NULL);
725 h->SetMarkerStyle(0);
726 h->SetBit(kCanDelete);
727 h->SetMarkerColor(kRed);
728 h->SetLineColor(kRed);
729 //h->SetFillColor(18);
730 h->Draw("same"/*"bar same"*/);
731
732 // This is the only way to make it work...
733 // Clone and copy constructor give strange crashes :(
734 h = fOffData->ProjectionZ("ProjOnOff", 0, fOffData->GetNbinsX()+1, 0, fOffData->GetNbinsY()+1, "E");
735 h->SetBit(TH1::kNoTitle);
736 h->SetXTitle(fHist.GetZaxis()->GetTitle());
737 h->SetYTitle("Counts");
738 h->SetDirectory(NULL);
739 h->SetMarkerStyle(0);
740 h->SetBit(kCanDelete);
741 h->SetMarkerColor(kBlue);
742 h->SetLineColor(kBlue);
743 h->Draw("same");
744
745 TLine lin;
746 lin.SetLineStyle(kDashed);
747 lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
748 }
749
750 // After the Alpha-projection has been drawn. Fit the histogram
751 // and paint the result into this pad
752 AppendPad("variable");
753
754 if (fHEnergy.GetNbinsX()>1 || fHEnergy.GetBinContent(1)>0)
755 {
756 pad->cd(2);
757 gPad->SetBorderMode(0);
758 gPad->SetGridx();
759 gPad->SetGridy();
760 fHEnergy.Draw();
761
762 AppendPad("energy");
763
764 h = (TH1D*)fHist.Project3D("y");
765 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
766 h->SetXTitle("E [GeV]");
767 h->SetYTitle("Counts");
768 h->SetDirectory(NULL);
769 h->SetMarkerStyle(kFullDotSmall);
770 h->SetBit(kCanDelete);
771 h->SetMarkerColor(kCyan);
772 h->SetLineColor(kCyan);
773 h->Draw("Psame");
774 }
775 else
776 delete pad->GetPad(2);
777
778 if ((fTimeEffOn && fTime) || fHTime.GetNbinsX()>1 || fHTime.GetBinError(1)>0)
779 {
780 pad->cd(3);
781 gPad->SetBorderMode(0);
782 gPad->SetGridx();
783 gPad->SetGridy();
784 fHTime.Draw();
785 AppendPad("time");
786 }
787 else
788 delete pad->GetPad(3);
789
790 if (fHTheta.GetNbinsX()>1 || fHTheta.GetBinContent(1)>0)
791 {
792 pad->cd(4);
793 gPad->SetGridx();
794 gPad->SetGridy();
795 gPad->SetBorderMode(0);
796 fHTheta.Draw();
797
798 AppendPad("theta");
799
800 h = (TH1D*)fHist.Project3D("x");
801 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
802 h->SetXTitle("\\theta [\\circ]");
803 h->SetYTitle("Counts");
804 h->SetDirectory(NULL);
805 h->SetMarkerStyle(kFullDotSmall);
806 h->SetBit(kCanDelete);
807 h->SetMarkerColor(kCyan);
808 h->SetLineColor(kCyan);
809 h->Draw("Psame");
810 }
811 else
812 delete pad->GetPad(4);
813}
814
815void MHAlpha::DrawAll(Bool_t newc)
816{
817 if (!newc && !fDisplay)
818 return;
819
820 // FIXME: Do in Paint if special option given!
821 TCanvas &c = newc ? *new TCanvas : fDisplay->AddTab("SizeBins");
822 Int_t n = fHist.GetNbinsY();
823 Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
824 c.Divide(nc, nc, 1e-10, 1e-10);
825 gPad->SetBorderMode(0);
826 gPad->SetFrameBorderMode(0);
827
828 // Do not store the 'final' result in fFit
829 MAlphaFitter fit(fFit);
830
831 for (int i=1; i<=fHist.GetNbinsY(); i++)
832 {
833 c.cd(i);
834 gPad->SetBorderMode(0);
835 gPad->SetFrameBorderMode(0);
836
837 TH1D *hon = fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, i, i, "E");
838 hon->SetBit(TH1::kNoTitle);
839 hon->SetStats(kTRUE);
840 hon->SetXTitle(fHist.GetZaxis()->GetTitle());
841 hon->SetYTitle("Counts");
842 hon->SetDirectory(NULL);
843 hon->SetMarkerStyle(0);
844 hon->SetBit(kCanDelete);
845 hon->Draw("");
846
847 TH1D *hof = 0;
848 Double_t alpha = 1;
849
850 if (fOffData)
851 {
852 hon->SetMarkerColor(kGreen);
853
854 hof = fOffData->ProjectionZ("ProjOff", 0, fOffData->GetNbinsX()+1, i, i, "E");
855 hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
856 hof->SetXTitle(fHist.GetZaxis()->GetTitle());
857 hof->SetYTitle("Counts");
858 hof->SetDirectory(NULL);
859 hof->SetMarkerStyle(0);
860 hof->SetBit(kCanDelete);
861 hof->SetMarkerColor(kRed);
862 hof->SetLineColor(kRed);
863 hof->Draw("same");
864
865 alpha = fit.Scale(*hof, *hon);
866
867 hon->SetMaximum();
868 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
869
870 TH1D *diff = new TH1D(*hon);
871 diff->Add(hof, -1);
872 diff->SetBit(TH1::kNoTitle);
873 diff->SetXTitle(fHist.GetZaxis()->GetTitle());
874 diff->SetYTitle("Counts");
875 diff->SetDirectory(NULL);
876 diff->SetMarkerStyle(0);
877 diff->SetBit(kCanDelete);
878 diff->SetMarkerColor(kBlue);
879 diff->SetLineColor(kBlue);
880 diff->Draw("same");
881
882 TLine lin;
883 lin.SetLineStyle(kDashed);
884 lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0);
885
886 const Float_t min = diff->GetMinimum()*1.05;
887 hon->SetMinimum(min<0 ? min : 0);
888 }
889
890 if (hof ? fit.Fit(*hon, *hof, alpha) : fit.Fit(*hon))
891 {
892 *fLog << dbg << "Bin " << i << ": sigmaexc=" << fit.GetEventsExcess()/fit.GetErrorExcess() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << " scale=" << fit.GetScaleFactor() << endl;
893 fit.DrawResult();
894 }
895 /*
896 if (fit.FitEnergy(fHist, fOffData, i, kTRUE))
897 {
898 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
899 fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
900 }*/
901 }
902
903}
904
905void MHAlpha::DrawNicePlot(Bool_t newc, const char *title, const char *watermark, Int_t binlo, Int_t binhi)
906{
907 if (!newc && !fDisplay)
908 return;
909
910 if (!fOffData)
911 return;
912
913 // Open and setup canvas/pad
914 TCanvas &c = newc ? *new TCanvas : fDisplay->AddTab("ThetsSq");
915
916 c.SetBorderMode(0);
917 c.SetFrameBorderMode(0);
918 c.SetFillColor(kWhite);
919
920 c.SetLeftMargin(0.12);
921 c.SetRightMargin(0.01);
922 c.SetBottomMargin(0.16);
923 c.SetTopMargin(0.18);
924
925 c.SetGridy();
926
927 gStyle->SetOptStat(0);
928
929 // Get on-data
930 TH1D *hon = (TH1D*)fHist.ProjectionZ("Proj", 0, fHist.GetNbinsX()+1, binlo, binhi, "E");
931 hon->SetDirectory(NULL);
932 hon->SetBit(kCanDelete);
933 hon->SetMarkerSize(0);
934 hon->SetLineWidth(2);
935 hon->SetLineColor(kBlack);
936 hon->SetMarkerColor(kBlack);
937
938 // Get off-data
939 TH1D *hoff = 0;
940 if (fOffData)
941 {
942 hoff = (TH1D*)fOffData->ProjectionZ("ProjOff", 0, fHist.GetNbinsX()+1, binlo, binhi, "E");
943 hoff->SetDirectory(NULL);
944 hoff->SetBit(kCanDelete);
945 hoff->SetFillColor(17);
946 hoff->SetMarkerSize(0);
947 hoff->SetLineColor(kBlack);
948 hoff->SetMarkerColor(kBlack);
949 }
950
951 // Setup plot which is drawn first
952 TH1D *h = hoff ? hoff : hon;
953 h->GetXaxis()->SetLabelSize(0.06);
954 h->GetXaxis()->SetTitleSize(0.06);
955 h->GetXaxis()->SetTitleOffset(0.95);
956 h->GetYaxis()->SetLabelSize(0.06);
957 h->GetYaxis()->SetTitleSize(0.06);
958 h->GetYaxis()->CenterTitle();
959 h->SetYTitle("Counts");
960 h->SetTitleSize(0.07);
961 h->SetTitle("");
962
963 const Double_t imax = fFit.GetSignalIntegralMax();
964// cout << "-----> IMAX " << imax << endl;
965// cout << " bins: " << binlo << " " << binhi << endl;
966// cout << " max: " << hon->GetMaximum() << " " << hoff->GetMaximum() << endl;
967// if (imax<1)
968// h->GetXaxis()->SetRangeUser(0, 0.6*0.6);
969
970 //to avoid that '0.3' is not displayed properly
971 h->GetXaxis()->SetRangeUser(0, 0.29);
972
973 // scale off-data
974
975 MAlphaFitter fit(fFit);
976 fit.ScaleAndFit(*hon, hoff);
977
978// cout << " max: " << hon->GetMaximum() << " " << hoff->GetMaximum() << endl;
979 hon->SetMinimum(0);
980 hoff->SetMinimum(0);
981
982 // draw data
983 if (hoff)
984 {
985 hoff->SetMaximum(TMath::Max(hon->GetMaximum(),hoff->GetMaximum())*1.1);
986// cout << "setting max to " << TMath::Max(hon->GetMaximum(),hoff->GetMaximum())*1.1 << endl;
987 hoff->Draw("bar");
988 hon->Draw("same");
989 }
990 else
991 {
992 hon->SetMaximum();
993 hon->Draw();
994 }
995
996 // draw a preliminary tag
997 TLatex text;
998 text.SetTextColor(kWhite);
999 text.SetBit(TLatex::kTextNDC);
1000 text.SetTextSize(0.07);
1001 text.SetTextAngle(2.5);
1002
1003 TString wm(watermark);
1004 if (binlo>=1 || binhi<hon->GetNbinsX())
1005 {
1006 wm += wm.IsNull() ? "(" : " (";
1007 if (binlo>=1)
1008 wm += MString::Format("%.1fGeV", fHist.GetYaxis()->GetBinLowEdge(binlo));
1009 wm += "-";
1010 if (binhi<hon->GetNbinsX())
1011 wm += MString::Format("%.1fGeV", fHist.GetYaxis()->GetBinLowEdge(binhi+1));
1012 wm += ")";
1013 }
1014 if (!wm.IsNull())
1015 text.DrawLatex(0.45, 0.2, wm);
1016 //enum { kTextNDC = BIT(14) };
1017
1018 // draw line showing cut
1019 TLine line;
1020 line.SetLineColor(14);
1021 line.SetLineStyle(7);
1022 line.DrawLine(imax, 0, imax, h->GetMaximum());
1023
1024 // Add a title above the plot
1025 TPaveText *pave=new TPaveText(0.12, 0.83, 0.99, 0.975, "blNDC");
1026 pave->SetBorderSize(1);
1027 pave->SetLabel(title);
1028
1029 TText *ptxt = pave->AddText(" ");
1030 ptxt->SetTextAlign(20);
1031
1032 ptxt = pave->AddText(MString::Format("Significance %.1f#sigma, off-scale %.2f",
1033 fit.GetSignificance(), fit.GetScaleFactor()));
1034// ptxt = pave->AddText(MString::Format("Significance %.1f, off-scale %.2f",
1035// fit.GetSignificance(), fit.GetScaleFactor()));
1036 ptxt->SetTextAlign(21);
1037 //ptxt->SetTextFont(4);
1038 //ptxt->SetTextSize(0.2);
1039
1040 ptxt = pave->AddText(MString::Format("%.1f excess events, %.1f background events",
1041 fit.GetEventsExcess(), fit.GetEventsBackground()));
1042
1043 ptxt->SetTextAlign(21);
1044 //ptxt->SetTextFont(4);
1045 //ptxt->SetTextSize(20);
1046
1047 pave->SetBit(kCanDelete);
1048 pave->Draw("br");
1049}
1050
1051Bool_t MHAlpha::Finalize()
1052{
1053 if (!FitAlpha())
1054 {
1055 *fLog << warn << "MAlphaFitter - Fit failed..." << endl;
1056 return kTRUE;
1057 }
1058
1059 // Store the final result in fFit
1060 *fLog << all;
1061 fFit.Print("result");
1062
1063 fResult->SetVal(fFit.GetMinimizationValue());
1064 fSigma->SetVal(fFit.GetGausSigma());
1065
1066 if (!fSkipHistEnergy)
1067 {
1068 *fLog << inf3 << "Processing energy bins..." << endl;
1069 FitEnergyBins();
1070 }
1071 if (!fSkipHistTheta)
1072 {
1073 *fLog << inf3 << "Processing theta bins..." << endl;
1074 FitThetaBins();
1075 }
1076 if (!fSkipHistTime)
1077 {
1078 *fLog << inf3 << "Processing time bins..." << endl;
1079 UpdateAlphaTime(kTRUE);
1080 MH::RemoveFirstBin(fHTime);
1081 }
1082
1083 if (fOffData)
1084 DrawAll(kFALSE);
1085
1086 return kTRUE;
1087}
1088
1089// --------------------------------------------------------------------------
1090//
1091// You can use this function if you want to use a MHMatrix instead of
1092// MMcEvt. This function adds all necessary columns to the
1093// given matrix. Afterward you should fill the matrix with the corresponding
1094// data (eg from a file by using MHMatrix::Fill). If you now loop
1095// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
1096// will take the values from the matrix instead of the containers.
1097//
1098// It takes fSkipHist* into account!
1099//
1100void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
1101{
1102 if (fMatrix)
1103 return;
1104
1105 fMatrix = mat;
1106
1107 fMap[0] = fMatrix->AddColumn(GetParameterRule());
1108 fMap[1] = -1;
1109 fMap[2] = -1;
1110 fMap[3] = -1;
1111 fMap[4] = -1;
1112
1113 if (!fSkipHistEnergy)
1114 {
1115 fMap[1] = type==0 ? fMatrix->AddColumn("MEnergyEst.fVal") : -1;
1116 fMap[2] = type==0 ? -1 : fMatrix->AddColumn("MHillas.fSize");
1117 }
1118
1119 if (!fSkipHistTheta)
1120 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
1121
1122 // if (!fSkipHistTime)
1123 // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
1124}
1125
1126void MHAlpha::StopMapping()
1127{
1128 fMatrix = NULL;
1129}
1130
1131void MHAlpha::ApplyScaling()
1132{
1133 if (!fOffData)
1134 return;
1135
1136 fFit.ApplyScaling(fHist, *const_cast<TH3D*>(fOffData));
1137}
1138
1139Int_t MHAlpha::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
1140{
1141 Bool_t rc = kFALSE;
1142 if (IsEnvDefined(env, prefix, "NumTimeBins", print))
1143 {
1144 SetNumTimeBins(GetEnvValue(env, prefix, "NumTimeBins", fNumTimeBins));
1145 rc = kTRUE;
1146 }
1147 if (IsEnvDefined(env, prefix, "ForceUsingSize", print))
1148 {
1149 fForceUsingSize = GetEnvValue(env, prefix, "ForceUsingSize", fForceUsingSize);
1150 rc = kTRUE;
1151 }
1152 return rc;
1153}
1154
1155Int_t MHAlpha::DistancetoPrimitive(Int_t px, Int_t py)
1156{
1157 // If pad has no subpad return (we are in one of the subpads)
1158 if (gPad->GetPad(1)==NULL)
1159 return 9999;
1160
1161 // If pad has a subpad we are in the main pad. Check its value.
1162 return gPad->GetPad(1)->DistancetoPrimitive(px,py)==0 ? 0 : 9999;
1163}
1164
1165void MHAlpha::RecursiveRemove(TObject *obj)
1166{
1167 if (obj==fOffData)
1168 fOffData = 0;
1169
1170 MH::RecursiveRemove(obj);
1171}
Note: See TracBrowser for help on using the repository browser.