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

Last change on this file since 8113 was 8016, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 27.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 <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 // For some reason the line-color is resetted
588 hoff->SetLineColor(kRed);
589
590 if ((h0=(TH1D*)gPad->FindObject("ProjOnOff")))
591 {
592 h0->Reset();
593 h0->Add(hoff, hon, -1);
594 const Float_t min = h0->GetMinimum()*1.05;
595 hon->SetMinimum(min<0 ? min : 0);
596 }
597 }
598 else
599 hon->SetMinimum(0);
600 FitEnergyBins();
601 FitThetaBins();
602 }
603
604 if (o==(TString)"variable")
605 if ((h0 = (TH1D*)gPad->FindObject("Proj")))
606 {
607 // Check whether we are dealing with on-off analysis
608 TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
609
610 // BE CARFEULL: This is a really weird workaround!
611 const Double_t scale = !hoff || hoff->GetMaximum()<0 ? 1 : hoff->GetMaximum();
612
613 // Do not store the 'final' result in fFit
614 MAlphaFitter fit(fFit);
615 fit.Fit(*h0, hoff, scale, kTRUE);
616 fit.PaintResult();
617 }
618
619 if (o==(TString)"time")
620 PaintText(fHTime);//.Integral(), 0);
621
622 if (o==(TString)"theta")
623 {
624 TH1 *h = (TH1*)gPad->FindObject(Form("%s_x", fHist.GetName()));
625 if (h)
626 {
627 TH1D *h2 = (TH1D*)fHist.Project3D("dum_x");
628 h2->SetDirectory(0);
629 h2->Scale(fHTheta.Integral()/h2->Integral());
630 h->Reset();
631 h->Add(h2);
632 delete h2;
633 }
634 PaintText(fHTheta);//.Integral(), 0);
635 }
636
637 if (o==(TString)"energy")
638 {
639 TH1 *h = (TH1*)gPad->FindObject(Form("%s_y", fHist.GetName()));
640 if (h)
641 {
642 TH1D *h2 = (TH1D*)fHist.Project3D("dum_y");
643 h2->SetDirectory(0);
644 h2->Scale(fHEnergy.Integral()/h2->Integral());
645 h->Reset();
646 h->Add(h2);
647 delete h2;
648 }
649 PaintText(fHEnergy);//.Integral(), 0);
650
651 if (fHEnergy.GetMaximum()>1)
652 {
653 gPad->SetLogx();
654 gPad->SetLogy();
655 }
656 }
657
658 gPad = padsave;
659}
660
661// --------------------------------------------------------------------------
662//
663// Draw the histogram
664//
665void MHAlpha::Draw(Option_t *opt)
666{
667 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
668
669 /*
670 if (TString(opt).Contains("sizebins", TString::kIgnoreCase))
671 {
672 AppendPad("sizebins");
673 return;
674 }
675 */
676
677 // Do the projection before painting the histograms into
678 // the individual pads
679 AppendPad("proj");
680
681 pad->SetBorderMode(0);
682 pad->Divide(2,2);
683
684 TH1D *h=0;
685
686 pad->cd(1);
687 gPad->SetBorderMode(0);
688
689 h = fHist.ProjectionZ("Proj", -1, 9999, -1, 9999, "E");
690 h->SetBit(TH1::kNoTitle);
691 h->SetStats(kTRUE);
692 h->SetXTitle(fHist.GetZaxis()->GetTitle());
693 h->SetYTitle("Counts");
694 h->SetDirectory(NULL);
695 h->SetMarkerStyle(kFullDotMedium);
696 h->SetBit(kCanDelete);
697 h->Draw("");
698
699 if (fOffData)
700 {
701 // To get green on-data
702 //h->SetMarkerColor(kGreen);
703 //h->SetLineColor(kGreen);
704
705 h = fOffData->ProjectionZ("ProjOff", -1, 9999, -1, 9999, "E");
706 h->SetBit(TH1::kNoTitle);
707 h->SetXTitle(fHist.GetZaxis()->GetTitle());
708 h->SetYTitle("Counts");
709 h->SetDirectory(NULL);
710 h->SetMarkerStyle(kFullDotMedium);
711 h->SetBit(kCanDelete);
712 h->SetMarkerColor(kRed);
713 h->SetLineColor(kRed);
714 //h->SetFillColor(18);
715 h->Draw("same"/*"bar same"*/);
716
717 h = (TH1D*)h->Clone("ProjOnOff");
718 h->SetBit(TH1::kNoTitle);
719 h->SetXTitle(fHist.GetZaxis()->GetTitle());
720 h->SetYTitle("Counts");
721 h->SetDirectory(NULL);
722 h->SetMarkerStyle(kFullDotMedium);
723 h->SetBit(kCanDelete);
724 h->SetMarkerColor(kBlue);
725 h->SetLineColor(kBlue);
726 h->Draw("same");
727
728 TLine lin;
729 lin.SetLineStyle(kDashed);
730 lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
731 }
732
733 // After the Alpha-projection has been drawn. Fit the histogram
734 // and paint the result into this pad
735 AppendPad("variable");
736
737 if (fHEnergy.GetNbinsX()>1)
738 {
739 pad->cd(2);
740 gPad->SetBorderMode(0);
741 gPad->SetGridx();
742 gPad->SetGridy();
743 fHEnergy.Draw();
744
745 AppendPad("energy");
746
747 h = (TH1D*)fHist.Project3D("y");
748 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
749 h->SetXTitle("E [GeV]");
750 h->SetYTitle("Counts");
751 h->SetDirectory(NULL);
752 h->SetMarkerStyle(kFullDotMedium);
753 h->SetBit(kCanDelete);
754 h->SetMarkerColor(kCyan);
755 h->SetLineColor(kCyan);
756 h->Draw("Psame");
757 }
758 else
759 delete pad->GetPad(2);
760
761 if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1)
762 {
763 pad->cd(3);
764 gPad->SetBorderMode(0);
765 gPad->SetGridx();
766 gPad->SetGridy();
767 fHTime.Draw();
768 AppendPad("time");
769 }
770 else
771 delete pad->GetPad(3);
772
773 if (fHTheta.GetNbinsX()>1)
774 {
775 pad->cd(4);
776 gPad->SetGridx();
777 gPad->SetGridy();
778 gPad->SetBorderMode(0);
779 fHTheta.Draw();
780
781 AppendPad("theta");
782
783 h = (TH1D*)fHist.Project3D("x");
784 h->SetBit(TH1::kNoTitle|TH1::kNoStats);
785 h->SetXTitle("\\theta [\\circ]");
786 h->SetYTitle("Counts");
787 h->SetDirectory(NULL);
788 h->SetMarkerStyle(kFullDotMedium);
789 h->SetBit(kCanDelete);
790 h->SetMarkerColor(kCyan);
791 h->SetLineColor(kCyan);
792 h->Draw("Psame");
793 }
794 else
795 delete pad->GetPad(4);
796}
797
798void MHAlpha::DrawAll(Bool_t newc)
799{
800 if (!newc && !fDisplay)
801 return;
802
803 // FIXME: Do in Paint if special option given!
804 TCanvas &c = newc ? *new TCanvas : fDisplay->AddTab("SizeBins");
805 Int_t n = fHist.GetNbinsY();
806 Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
807 c.Divide(nc, nc, 1e-10, 1e-10);
808 gPad->SetBorderMode(0);
809 gPad->SetFrameBorderMode(0);
810
811 // Do not store the 'final' result in fFit
812 MAlphaFitter fit(fFit);
813
814 for (int i=1; i<=fHist.GetNbinsY(); i++)
815 {
816 c.cd(i);
817 gPad->SetBorderMode(0);
818 gPad->SetFrameBorderMode(0);
819
820 TH1D *hon = fHist.ProjectionZ("Proj", -1, 9999, i, i, "E");
821 hon->SetBit(TH1::kNoTitle);
822 hon->SetStats(kTRUE);
823 hon->SetXTitle(fHist.GetZaxis()->GetTitle());
824 hon->SetYTitle("Counts");
825 hon->SetDirectory(NULL);
826 hon->SetMarkerStyle(kFullDotMedium);
827 hon->SetBit(kCanDelete);
828 hon->Draw("");
829
830 TH1D *hof = 0;
831 Double_t alpha = 1;
832
833 if (fOffData)
834 {
835 hon->SetMarkerColor(kGreen);
836
837 hof = fOffData->ProjectionZ("ProjOff", -1, 9999, i, i, "E");
838 hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
839 hof->SetXTitle(fHist.GetZaxis()->GetTitle());
840 hof->SetYTitle("Counts");
841 hof->SetDirectory(NULL);
842 hof->SetMarkerStyle(kFullDotMedium);
843 hof->SetBit(kCanDelete);
844 hof->SetMarkerColor(kRed);
845 hof->SetLineColor(kRed);
846 hof->Draw("same");
847
848 alpha = fit.Scale(*hof, *hon);
849
850 hon->SetMaximum();
851 hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
852
853 TH1D *diff = new TH1D(*hon);
854 diff->Add(hof, -1);
855 diff->SetBit(TH1::kNoTitle);
856 diff->SetXTitle(fHist.GetZaxis()->GetTitle());
857 diff->SetYTitle("Counts");
858 diff->SetDirectory(NULL);
859 diff->SetMarkerStyle(kFullDotMedium);
860 diff->SetBit(kCanDelete);
861 diff->SetMarkerColor(kBlue);
862 diff->Draw("same");
863
864 TLine lin;
865 lin.SetLineStyle(kDashed);
866 lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0);
867
868 const Float_t min = diff->GetMinimum()*1.05;
869 hon->SetMinimum(min<0 ? min : 0);
870 }
871
872 if (hof ? fit.Fit(*hon, *hof, alpha) : fit.Fit(*hon))
873 {
874 *fLog << dbg << "Bin " << i << ": sigmaexc=" << fit.GetSignificanceExc() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << " scale=" << fit.GetScaleFactor() << endl;
875 fit.DrawResult();
876 }
877 /*
878 if (fit.FitEnergy(fHist, fOffData, i, kTRUE))
879 {
880 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
881 fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
882 }*/
883 }
884
885}
886
887
888Bool_t MHAlpha::Finalize()
889{
890 //TH1D *h = fHist.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
891 //h->SetDirectory(0);
892 //Bool_t rc = fFit.Fit(*h);
893 //delete h;
894
895 if (!fFit.FitAlpha(fHist, fOffData))
896 {
897 *fLog << warn << "MAlphaFitter - Fit failed..." << endl;
898 return kTRUE;
899 }
900
901 // Store the final result in fFit
902 *fLog << all;
903 fFit.Print("result");
904
905 fResult->SetVal(fFit.GetMinimizationValue());
906 fSigma->SetVal(fFit.GetGausSigma());
907
908 if (!fSkipHistEnergy)
909 {
910 *fLog << inf << "Processing energy bins..." << endl;
911 FitEnergyBins();
912 }
913 if (!fSkipHistTheta)
914 {
915 *fLog << inf << "Processing theta bins..." << endl;
916 FitThetaBins();
917 }
918 if (!fSkipHistTime)
919 {
920 *fLog << inf << "Processing time bins..." << endl;
921 UpdateAlphaTime(kTRUE);
922 MH::RemoveFirstBin(fHTime);
923 }
924
925 if (fOffData)
926 DrawAll(kFALSE);
927
928 return kTRUE;
929}
930
931// --------------------------------------------------------------------------
932//
933// You can use this function if you want to use a MHMatrix instead of
934// MMcEvt. This function adds all necessary columns to the
935// given matrix. Afterward you should fill the matrix with the corresponding
936// data (eg from a file by using MHMatrix::Fill). If you now loop
937// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
938// will take the values from the matrix instead of the containers.
939//
940// It takes fSkipHist* into account!
941//
942void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
943{
944 if (fMatrix)
945 return;
946
947 fMatrix = mat;
948
949 fMap[0] = fMatrix->AddColumn(GetParameterRule());
950 fMap[1] = -1;
951 fMap[2] = -1;
952 fMap[3] = -1;
953 fMap[4] = -1;
954
955 if (!fSkipHistEnergy)
956 if (type==0)
957 {
958 fMap[1] = fMatrix->AddColumn("MEnergyEst.fVal");
959 fMap[2] = -1;
960 }
961 else
962 {
963 fMap[1] = -1;
964 fMap[2] = fMatrix->AddColumn("MHillas.fSize");
965 }
966
967 if (!fSkipHistTheta)
968 fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
969
970 // if (!fSkipHistTime)
971 // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
972}
973
974void MHAlpha::StopMapping()
975{
976 fMatrix = NULL;
977}
978
979Int_t MHAlpha::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
980{
981 Bool_t rc = kFALSE;
982 if (IsEnvDefined(env, prefix, "NumTimeBins", print))
983 {
984 SetNumTimeBins(GetEnvValue(env, prefix, "NumTimeBins", fNumTimeBins));
985 rc = kTRUE;
986 }
987 if (IsEnvDefined(env, prefix, "ForceUsingSize", print))
988 {
989 fForceUsingSize = GetEnvValue(env, prefix, "ForceUsingSize", fForceUsingSize);
990 rc = kTRUE;
991 }
992 return rc;
993}
994
995Int_t MHAlpha::DistancetoPrimitive(Int_t px, Int_t py)
996{
997 // If pad has no subpad return (we are in one of the subpads)
998 if (gPad->GetPad(1)==NULL)
999 return 9999;
1000
1001 // If pad has a subpad we are in the main pad. Check its value.
1002 return gPad->GetPad(1)->DistancetoPrimitive(px,py)==0 ? 0 : 9999;
1003}
Note: See TracBrowser for help on using the repository browser.