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

Last change on this file since 5098 was 5012, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 14.7 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-2004
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#include "MHAlpha.h"
40
41#include <TF1.h>
42#include <TGraph.h>
43#include <TStyle.h>
44#include <TLatex.h>
45#include <TCanvas.h>
46#include <TPaveStats.h>
47#include <TStopwatch.h>
48
49#include "MGeomCam.h"
50#include "MSrcPosCam.h"
51#include "MHillasSrc.h"
52#include "MEnergyEst.h"
53#include "MTime.h"
54#include "MObservatory.h"
55#include "MPointingPos.h"
56#include "MAstroSky2Local.h"
57#include "MStatusDisplay.h"
58#include "MParameters.h"
59#include "MHMatrix.h"
60
61#include "MBinning.h"
62#include "MParList.h"
63
64#include "MLog.h"
65#include "MLogManip.h"
66
67ClassImp(MHAlpha);
68ClassImp(MAlphaFitter);
69
70using namespace std;
71
72// --------------------------------------------------------------------------
73//
74// Default Constructor
75//
76MHAlpha::MHAlpha(const char *name, const char *title)
77 : fResult(0), /*fExcess(0),*/ fEnergy(0), fPointPos(0), fTimeEffOn(0),
78 fTime(0), fNameProjAlpha(Form("Alpha_%p", this)), fMatrix(0)
79{
80 //
81 // set the name and title of this object
82 //
83 fName = name ? name : "MHAlpha";
84 fTitle = title ? title : "Alpha plot";
85
86 fHAlpha.SetName("Alpha");
87 fHAlpha.SetTitle("Alpha");
88 fHAlpha.SetXTitle("\\Theta [deg]");
89 fHAlpha.SetYTitle("E_{est} [GeV]");
90 fHAlpha.SetZTitle("|\\alpha| [\\circ]");
91 fHAlpha.SetDirectory(NULL);
92 fHAlpha.UseCurrentStyle();
93
94 // Main histogram
95 fHAlphaTime.SetName("Alpha");
96 fHAlphaTime.SetXTitle("|\\alpha| [\\circ]");
97 fHAlphaTime.SetYTitle("Counts");
98 fHAlphaTime.UseCurrentStyle();
99 fHAlphaTime.SetDirectory(NULL);
100
101
102 fHEnergy.SetName("Energy");
103 fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
104 fHEnergy.SetXTitle("E_{est} [GeV]");
105 fHEnergy.SetYTitle("N_{exc}");
106 fHEnergy.SetDirectory(NULL);
107 fHEnergy.UseCurrentStyle();
108
109 fHTheta.SetName("Theta");
110 fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
111 fHTheta.SetXTitle("\\Theta [\\circ]");
112 fHTheta.SetYTitle("N_{exc}");
113 fHTheta.SetDirectory(NULL);
114 fHTheta.UseCurrentStyle();
115
116 // effective on time versus time
117 fHTime.SetName("Time");
118 fHTime.SetTitle(" N_{exc} vs. Time ");
119 fHTime.SetXTitle("Time");
120 fHTime.SetYTitle("N_{exc} [s]");
121 fHTime.UseCurrentStyle();
122 fHTime.SetDirectory(NULL);
123 fHTime.GetYaxis()->SetTitleOffset(1.2);
124 fHTime.GetXaxis()->SetLabelSize(0.033);
125 fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
126 fHTime.GetXaxis()->SetTimeDisplay(1);
127 fHTime.Sumw2();
128
129 MBinning binsa, binse, binst;
130 binsa.SetEdges(18, 0, 90);
131 binse.SetEdgesLog(25, 10, 100000);
132 binst.SetEdgesCos(50, 0, 60);
133 binse.Apply(fHEnergy);
134 binst.Apply(fHTheta);
135 binsa.Apply(fHAlphaTime);
136
137 MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
138}
139
140void MHAlpha::FitEnergySpec(Bool_t paint)
141{
142 /*
143 TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]*x");//
144 f.SetParameter(0, -2.1);
145 f.SetParameter(1, 1400); // 50
146 f.SetParameter(2, fHEnergy.Integral()); // 50
147 f.SetLineWidth(2);
148
149 fHEnergy.Fit(&f, "0QI", "", 100, 2400); // Integral? 90, 2800
150
151 const Double_t chi2 = f.GetChisquare();
152 const Int_t NDF = f.GetNDF();
153
154 // was fit successful ?
155 const Bool_t ok = NDF>0 && chi2<2.5*NDF;
156 f.SetLineColor(ok ? kGreen : kRed);
157
158 if (paint)
159 {
160 f.Paint("same");
161
162 if (ok)
163 {
164 TString s;
165 s += Form("\\alpha=%.1f ", f.GetParameter(0));
166 s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000);
167 s += Form("N=%.1e", f.GetParameter(2));
168
169 TLatex text(0.25, 0.94, s.Data());
170 text.SetBit(TLatex::kTextNDC);
171 text.SetTextSize(0.04);
172 text.Paint();
173 }
174 }*/
175}
176
177void MHAlpha::FitEnergyBins(Bool_t paint)
178{
179 // Do not store the 'final' result in fFit
180 MAlphaFitter fit(fFit);
181
182 const Int_t n = fHAlpha.GetNbinsY();
183
184 for (int i=1; i<=n; i++)
185 {
186 TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", -1, 9999, i, i, i==1?"E":"");
187 if (fit.Fit(*h))
188 {
189 fHEnergy.SetBinContent(i, fit.GetEventsExcess());
190 fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
191 }
192 delete h;
193 }
194}
195
196void MHAlpha::FitThetaBins(Bool_t paint)
197{
198 // Do not store the 'final' result in fFit
199 MAlphaFitter fit(fFit);
200
201 const Int_t n = fHAlpha.GetNbinsX();
202
203 for (int i=1; i<=n; i++)
204 {
205 TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", i, i, -1, 9999, i==1?"E":"");
206 if (fit.Fit(*h))
207 {
208 fHTheta.SetBinContent(i, fit.GetEventsExcess());
209 fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2);
210 }
211 delete h;
212 }
213}
214
215Bool_t MHAlpha::SetupFill(const MParList *pl)
216{
217 fHAlpha.Reset();
218 fHEnergy.Reset();
219 fHTheta.Reset();
220
221 fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst");
222 if (!fEnergy)
223 *fLog << warn << "MEnergyEst not found... ignored." << endl;
224
225 fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
226 if (!fPointPos)
227 *fLog << warn << "MPointingPos not found... ignored." << endl;
228
229 fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
230 if (!fTimeEffOn)
231 *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
232
233 fTime = (MTime*)pl->FindObject("MTime");
234 if (!fTime)
235 *fLog << warn << "MTime not found... ignored." << endl;
236
237 fResult = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "Significance");
238 if (!fResult)
239 return kFALSE;
240
241 //fExcess = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MExcess");
242 //if (!fExcess)
243 // return kFALSE;
244
245 fLastTime = MTime();
246
247 MBinning binst, binse, binsa;
248 binst.SetEdges(fHAlpha, 'x');
249 binse.SetEdges(fHAlpha, 'y');
250 binsa.SetEdges(fHAlpha, 'z');
251
252 MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning");
253 if (fPointPos && bins)
254 binst.SetEdges(*bins);
255 if (!fPointPos)
256 binst.SetEdges(1, 0, 90);
257
258 bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
259 if (fEnergy && bins)
260 binse.SetEdges(*bins);
261 if (!fEnergy)
262 binse.SetEdges(1, 10, 100000);
263
264 bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning");
265 if (bins)
266 binsa.SetEdges(*bins);
267
268 binse.Apply(fHEnergy);
269 binst.Apply(fHTheta);
270 binsa.Apply(fHAlphaTime);
271 MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
272
273 MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
274 if (!fit)
275 *fLog << warn << "MAlphaFitter not found... ignored." << endl;
276 else
277 fit->Copy(fFit);
278
279 return kTRUE;
280}
281
282void MHAlpha::InitAlphaTime(const MTime &t)
283{
284 //
285 // If this is the first call we have to initialize the time-histogram
286 //
287 MBinning bins;
288 bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime());
289 bins.Apply(fHTime);
290
291 fLastTime=t;
292}
293
294void MHAlpha::UpdateAlphaTime(Bool_t final)
295{
296 if (!fTimeEffOn)
297 return;
298
299 const Int_t steps = 6;
300
301 static int rebin = steps;
302
303 if (!final)
304 {
305 if (fLastTime==MTime() && *fTimeEffOn!=MTime())
306 {
307 InitAlphaTime(*fTimeEffOn);
308 return;
309 }
310
311 if (fLastTime!=*fTimeEffOn)
312 {
313 fLastTime=*fTimeEffOn;
314 rebin--;
315 }
316
317 if (rebin>0)
318 return;
319 }
320
321 MAlphaFitter fit(fFit);
322 if (!fit.Fit(fHAlphaTime))
323 return;
324
325 // Reset Histogram
326 fHAlphaTime.Reset();
327
328 // Get number of bins
329 const Int_t n = fHTime.GetNbinsX();
330
331 //
332 // Prepare Histogram
333 //
334
335 MTime *time = final ? fTime : fTimeEffOn;
336 if (final)
337 time->Plus1ns();
338
339 // Enhance binning
340 MBinning bins;
341 bins.SetEdges(fHTime, 'x');
342 bins.AddEdge(time->GetAxisTime());
343 bins.Apply(fHTime);
344
345 //
346 // Fill histogram
347 //
348 fHTime.SetBinContent(n+1, fit.GetEventsExcess());
349 fHTime.SetBinError(n+1, fit.GetEventsExcess()*0.1);
350
351 *fLog << all << *fTimeEffOn << ": " << fit.GetEventsExcess() << endl;
352
353 rebin = steps;
354}
355
356// --------------------------------------------------------------------------
357//
358// Fill the histogram. For details see the code or the class description
359//
360Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
361{
362 Double_t alpha, energy, theta;
363
364 if (fMatrix)
365 {
366 alpha = (*fMatrix)[fMap[0]];
367 energy = 1000;
368 theta = 0;
369 }
370 else
371 {
372 const MHillasSrc *hil = dynamic_cast<const MHillasSrc*>(par);
373 if (!par)
374 {
375 *fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
376 return kFALSE;
377 }
378
379 alpha = hil->GetAlpha();
380 energy = fEnergy ? fEnergy->GetEnergy() : 1000;
381 theta = fPointPos ? fPointPos->GetZd() : 0;
382 }
383
384 UpdateAlphaTime();
385
386 fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
387 fHAlphaTime.Fill(TMath::Abs(alpha), w);
388
389 return kTRUE;
390}
391
392// --------------------------------------------------------------------------
393//
394// Paint the integral and the error on top of the histogram
395//
396void MHAlpha::PaintText(Double_t val, Double_t error) const
397{
398 TLatex text(0.45, 0.94, Form("N_{exc} = %.1fs \\pm %.1fs", val, error));
399 text.SetBit(TLatex::kTextNDC);
400 text.SetTextSize(0.04);
401 text.Paint();
402}
403
404// --------------------------------------------------------------------------
405//
406// Update the projections
407//
408void MHAlpha::Paint(Option_t *opt)
409{
410 // Note: Any object cannot be added to a pad twice!
411 // The object is searched by gROOT->FindObject only in
412 // gPad itself!
413 TVirtualPad *padsave = gPad;
414
415 TH1D *h0=0;
416
417 TString o(opt);
418 if (o==(TString)"proj")
419 {
420 TPaveStats *st=0;
421 for (int x=0; x<4; x++)
422 {
423 TVirtualPad *p=padsave->GetPad(x+1);
424 if (!p || !(st = (TPaveStats*)p->GetPrimitive("stats")))
425 continue;
426
427 if (st->GetOptStat()==11)
428 continue;
429
430 const Double_t y1 = st->GetY1NDC();
431 const Double_t y2 = st->GetY2NDC();
432 const Double_t x1 = st->GetX1NDC();
433 const Double_t x2 = st->GetX2NDC();
434
435 st->SetY1NDC((y2-y1)/3+y1);
436 st->SetX1NDC((x2-x1)/3+x1);
437 st->SetOptStat(11);
438 }
439
440 padsave->cd(1);
441 fHAlpha.ProjectionZ(fNameProjAlpha);
442 FitEnergyBins();
443 FitThetaBins();
444 }
445
446 if (o==(TString)"alpha")
447 if ((h0 = (TH1D*)gPad->FindObject(fNameProjAlpha)))
448 {
449 // Do not store the 'final' result in fFit
450 MAlphaFitter fit(fFit);
451 fit.Fit(*h0, kTRUE);
452 fit.PaintResult();
453 }
454
455 if (o==(TString)"time")
456 PaintText(fHTime.Integral(), 0);
457
458 if (o==(TString)"theta")
459 PaintText(fHTheta.Integral(), 0);
460
461 if (o==(TString)"energy")
462 {
463 if (fHEnergy.GetEntries()>10)
464 {
465 gPad->SetLogx();
466 gPad->SetLogy();
467 }
468 FitEnergySpec(kTRUE);
469
470 }
471
472 gPad = padsave;
473}
474
475// --------------------------------------------------------------------------
476//
477// Draw the histogram
478//
479void MHAlpha::Draw(Option_t *opt)
480{
481 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
482
483 // Do the projection before painting the histograms into
484 // the individual pads
485 AppendPad("proj");
486
487 pad->SetBorderMode(0);
488 pad->Divide(2,2);
489
490 TH1D *h=0;
491
492 pad->cd(1);
493 gPad->SetBorderMode(0);
494 h = fHAlpha.ProjectionZ(fNameProjAlpha, -1, 9999, -1, 9999, "E");
495 h->SetBit(TH1::kNoTitle);
496 h->SetXTitle("\\alpha [\\circ]");
497 h->SetYTitle("Counts");
498 h->SetDirectory(NULL);
499 h->SetMarkerStyle(kFullDotMedium);
500 h->SetBit(kCanDelete);
501 h->Draw();
502 // After the Alpha-projection has been drawn. Fit the histogram
503 // and paint the result into this pad
504 AppendPad("alpha");
505
506 if (fHEnergy.GetNbinsX()>1)
507 {
508 pad->cd(2);
509 gPad->SetBorderMode(0);
510 fHEnergy.Draw();
511 AppendPad("energy");
512 }
513 else
514 delete pad->GetPad(2);
515
516 if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1)
517 {
518 pad->cd(3);
519 gPad->SetBorderMode(0);
520 fHTime.Draw();
521 AppendPad("time");
522 }
523 else
524 delete pad->GetPad(3);
525
526 if (fHTheta.GetNbinsX()>1)
527 {
528 pad->cd(4);
529 gPad->SetBorderMode(0);
530 fHTheta.Draw();
531 AppendPad("theta");
532 }
533 else
534 delete pad->GetPad(4);
535
536}
537
538Bool_t MHAlpha::Finalize()
539{
540 // Store the final result in fFit
541 TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
542 Bool_t rc = fFit.Fit(*h);
543 delete h;
544 if (!rc)
545 {
546 *fLog << warn << "Histogram empty." << endl;
547 return kTRUE;
548 }
549
550 if (fResult)
551 fResult->SetVal(fFit.GetSignificance());
552
553 FitEnergyBins();
554 FitThetaBins();
555 UpdateAlphaTime(kTRUE);
556 MH::RemoveFirstBin(fHTime);
557
558 return kTRUE;
559}
560
561// --------------------------------------------------------------------------
562//
563// You can use this function if you want to use a MHMatrix instead of
564// MMcEvt. This function adds all necessary columns to the
565// given matrix. Afterward you should fill the matrix with the corresponding
566// data (eg from a file by using MHMatrix::Fill). If you now loop
567// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
568// will take the values from the matrix instead of the containers.
569//
570void MHAlpha::InitMapping(MHMatrix *mat)
571{
572 if (fMatrix)
573 return;
574
575 fMatrix = mat;
576
577 fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
578 //fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
579}
580
581void MHAlpha::StopMapping()
582{
583 fMatrix = NULL;
584}
Note: See TracBrowser for help on using the repository browser.