source: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc@ 7825

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