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

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