source: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc@ 7169

Last change on this file since 7169 was 7169, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 29.6 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, 4/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJSpectrum
28//
29// Program to calculate spectrum
30//
31/////////////////////////////////////////////////////////////////////////////
32#include "MJSpectrum.h"
33
34// Root
35#include <TF1.h>
36#include <TH1.h>
37#include <TH2.h>
38#include <TFile.h>
39#include <TChain.h>
40#include <TLatex.h>
41#include <TCanvas.h>
42#include <TObjArray.h>
43
44// Environment
45#include "MLog.h"
46#include "MLogManip.h"
47
48#include "MStatusArray.h"
49#include "MStatusDisplay.h"
50
51// Container
52#include "MH3.h"
53#include "MBinning.h"
54#include "MDataSet.h"
55#include "MMcCorsikaRunHeader.h"
56
57// Spectrum
58#include "../mhflux/MAlphaFitter.h"
59#include "../mhflux/MHAlpha.h"
60#include "../mhflux/MHCollectionArea.h"
61#include "../mhflux/MHEnergyEst.h"
62#include "../mhflux/MMcSpectrumWeight.h"
63
64// Eventloop
65#include "MEvtLoop.h"
66#include "MTaskList.h"
67#include "MParList.h"
68
69// Tasks/Filter
70#include "MReadMarsFile.h"
71#include "MReadMarsFile.h"
72#include "MFEventSelector2.h"
73#include "MFDataMember.h"
74#include "MEnergyEstimate.h"
75#include "MTaskEnv.h"
76#include "MFillH.h"
77#include "MHillasCalc.h"
78#include "MSrcPosCalc.h"
79#include "MContinue.h"
80
81ClassImp(MJSpectrum);
82
83using namespace std;
84
85MJSpectrum::MJSpectrum(const char *name, const char *title)
86 : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
87 fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE),
88 fNoThetaWeights(kFALSE)
89{
90 fName = name ? name : "MJSpectrum";
91 fTitle = title ? title : "Standard program to calculate spectrum";
92}
93
94MJSpectrum::~MJSpectrum()
95{
96 if (fCut0)
97 delete fCut0;
98 if (fCut1)
99 delete fCut1;
100 if (fCut2)
101 delete fCut2;
102 if (fCut3)
103 delete fCut3;
104 if (fEstimateEnergy)
105 delete fEstimateEnergy;
106}
107
108// --------------------------------------------------------------------------
109//
110// Setup a task estimating the energy. The given task is cloned.
111//
112void MJSpectrum::SetEnergyEstimator(const MTask *task)
113{
114 if (fEstimateEnergy)
115 delete fEstimateEnergy;
116 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
117}
118
119Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const
120{
121 if (task)
122 {
123 delete task;
124 task = 0;
125 }
126
127 task = (MTask*)gFile->Get(name);
128 if (!task)
129 {
130 *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
131 return kFALSE;
132 }
133 if (!task->InheritsFrom(MTask::Class()))
134 {
135 *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
136 delete task;
137 return kFALSE;
138 }
139
140 task->SetName(name);
141
142 return kTRUE;
143}
144
145void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
146{
147 fLog->Separator("Alpha Fitter");
148 *fLog << all;
149 fit.Print();
150
151 fLog->Separator("Used Cuts");
152 fCut0->Print();
153 fCut1->Print();
154 fCut2->Print();
155 fCut3->Print();
156
157 //gLog.Separator("Energy Estimator");
158 //fEstimateEnergy->Print();
159}
160
161// --------------------------------------------------------------------------
162//
163// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
164// the dataset.
165// The simulated energy range and spectral slope is initialized from
166// there.
167// In the following eventloops the forced check in MMcSpectrumWeight
168// ensures, that the spectral slope and energy range doesn't change.
169//
170Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
171{
172 fLog->Separator("Initialize energy weighting");
173
174 if (!CheckEnv(w))
175 {
176 *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
177 return kFALSE;
178 }
179
180 TChain chain("RunHeaders");
181 set.AddFilesOn(chain);
182
183 MMcCorsikaRunHeader *h=0;
184 chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
185 chain.GetEntry(1);
186
187 if (!h)
188 {
189 *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
190 return kFALSE;
191 }
192
193 if (!w.Set(*h))
194 {
195 *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
196 return kFALSE;
197 }
198
199 w.Print();
200 return kTRUE;
201}
202
203Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
204{
205 *fLog << inf << "Reading from file: " << fPathIn << endl;
206
207 TFile file(fPathIn, "READ");
208 if (!file.IsOpen())
209 {
210 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
211 return -1;
212 }
213
214 MStatusArray arr;
215 if (arr.Read()<=0)
216 {
217 *fLog << "MStatusDisplay not found in file... abort." << endl;
218 return -1;
219 }
220
221 TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
222 TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
223 if (!vstime || !size)
224 return -1;
225
226 vstime->Copy(h1);
227 size->Copy(h2);
228 h1.SetDirectory(0);
229 h2.SetDirectory(0);
230
231 if (fDisplay)
232 arr.DisplayIn(*fDisplay, "Hist");
233
234 if (!ReadTask(fCut0, "Cut0"))
235 return -1;
236 if (!ReadTask(fCut1, "Cut1"))
237 return -1;
238 if (!ReadTask(fCut2, "Cut2"))
239 return -1;
240 if (!ReadTask(fCut3, "Cut3"))
241 return -1;
242
243 TObjArray arrread;
244
245 TIter Next(plist);
246 TObject *o=0;
247 while ((o=Next()))
248 if (o->InheritsFrom(MBinning::Class()))
249 arrread.Add(o);
250
251 arrread.Add(plist.FindObject("MAlphaFitter"));
252
253 if (!ReadContainer(arrread))
254 return -1;
255
256 return vstime->Integral();
257}
258
259Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
260{
261 // Some debug output
262 fLog->Separator("Compiling original MC distribution");
263
264 weight.SetNameMcEvt("MMcEvtBasic");
265 const TString w(weight.GetFormulaWeights());
266 weight.SetNameMcEvt();
267
268 *fLog << inf << "Using weights: " << w << endl;
269 *fLog << "Please stand by, this may take a while..." << flush;
270
271 if (fDisplay)
272 fDisplay->SetStatusLine1("Compiling MC distribution...");
273
274 // Create chain
275 TChain chain("OriginalMC");
276 set.AddFilesOn(chain);
277
278 // Prepare histogram
279 h.Reset();
280
281 // Fill histogram from chain
282 h.SetDirectory(gROOT);
283 if (h.InheritsFrom(TH2::Class()))
284 {
285 h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
286 h.SetXTitle("\\Theta [\\circ]");
287 h.SetYTitle("E [GeV]");
288 h.SetZTitle("Counts");
289 chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
290 }
291 else
292 {
293 h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
294 h.SetXTitle("\\Theta [\\circ]");
295 h.SetYTitle("Counts");
296 chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
297 }
298 h.SetDirectory(0);
299
300 *fLog << "done." << endl;
301 if (fDisplay)
302 fDisplay->SetStatusLine2("done.");
303
304 if (h.GetEntries()>0)
305 return kTRUE;
306
307 *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
308
309 return h.GetEntries()>0;
310}
311
312Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
313{
314 // Display some stuff
315 if (fDisplay)
316 {
317 TCanvas &c = fDisplay->AddTab("ZdDist");
318 c.Divide(2,2);
319
320 // On-Time vs. Theta
321 c.cd(1);
322 gPad->SetBorderMode(0);
323 temp1.DrawCopy();
324
325 // Number of MC events (produced) vs Theta
326 c.cd(2);
327 gPad->SetBorderMode(0);
328 temp2.SetName("NVsTheta");
329 temp2.DrawCopy();
330
331 c.cd(4);
332 gPad->SetBorderMode(0);
333
334 c.cd(3);
335 gPad->SetBorderMode(0);
336 }
337
338 // Calculate the Probability
339 temp1.Divide(&temp2);
340 temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral());
341
342 // Some cosmetics: Name, Axis, etc.
343 temp1.SetName("ProbVsTheta");
344 temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
345 temp1.SetYTitle("Probability");
346 if (fDisplay)
347 temp1.DrawCopy();
348
349 return kTRUE;
350}
351
352// --------------------------------------------------------------------------
353//
354// Display the final theta distribution.
355//
356void MJSpectrum::DisplayResult(const TH2D &h2) const
357{
358 if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
359 return;
360
361 TH1D &proj = *h2.ProjectionX();
362 proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
363 proj.SetXTitle("\\Theta [\\circ]");
364 proj.SetYTitle("Counts");
365 proj.SetLineColor(kBlue);
366 proj.SetDirectory(0);
367 proj.SetBit(kCanDelete);
368
369 TVirtualPad *pad = gPad;
370
371 pad->cd(4);
372 proj.DrawCopy();
373
374 pad->cd(1);
375 TH1D *theta = (TH1D*)gPad->FindObject("Theta");
376 if (theta)
377 {
378 proj.Scale(theta->GetMaximum()/proj.GetMaximum());
379 theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
380 }
381 proj.Draw("same");
382}
383
384// --------------------------------------------------------------------------
385//
386// Fills the excess histogram (vs E-est) from the events stored in the
387// ganymed result file and therefor estimates the energy.
388//
389// The resulting histogram excess-vs-energy ist copied into h2.
390//
391Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
392{
393 // Try to find the class used to determin the signal!
394 TString cls("MHAlpha");
395 if (fDisplay)
396 {
397 TCanvas *c = fDisplay->GetCanvas("Hist");
398 if (c)
399 {
400 TIter Next(c->GetListOfPrimitives());
401 TObject *obj=0;
402 while ((obj=Next()))
403 if (obj->InheritsFrom(MHAlpha::Class()))
404 break;
405 if (obj)
406 cls = obj->ClassName();
407 }
408 }
409
410 cout << "FOUND: "<< cls << endl;
411
412 // Now fill the histogram
413 *fLog << endl;
414 fLog->Separator("Refill Excess");
415 *fLog << endl;
416
417 MTaskList tlist;
418 plist.AddToList(&tlist);
419
420 MReadTree read("Events");
421 read.DisableAutoScheme();
422 read.AddFile(fPathIn);
423
424 MEnergyEstimate est;
425 MTaskEnv taskenv1("EstimateEnergy");
426 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
427
428 MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
429 MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
430
431 MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
432 MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
433
434 fill1.SetFilter(&f0);
435 fill2.SetFilter(&f1);
436
437 tlist.AddToList(&read);
438 tlist.AddToList(&taskenv1);
439 tlist.AddToList(&f0);
440 tlist.AddToList(&f1);
441 tlist.AddToList(&fill1);
442 tlist.AddToList(&fill2);
443
444 MEvtLoop loop(fName);
445 loop.SetParList(&plist);
446 loop.SetDisplay(fDisplay);
447 loop.SetLogStream(fLog);
448
449 if (!SetupEnv(loop))
450 return kFALSE;
451
452 if (!loop.Eventloop())
453 {
454 *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
455 return kFALSE;
456 }
457
458 if (!loop.GetDisplay())
459 {
460 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
461 return kFALSE;
462 }
463
464 const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE");
465 if (!halpha)
466 {
467 *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
468 return kFALSE;
469 }
470
471 halpha->GetHEnergy().Copy(h2);
472 h2.SetDirectory(0);
473
474 return kTRUE;
475}
476
477Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
478{
479 MTaskList tlist1;
480 plist.Replace(&tlist1);
481
482 MReadMarsFile readmc("OriginalMC");
483 //readmc.DisableAutoScheme();
484 set.AddFilesOn(readmc);
485 readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
486 readmc.EnableBranch("MMcEvtBasic.fEnergy");
487
488 mh1.SetLogy();
489 mh1.SetLogz();
490 mh1.SetName("ThetaE");
491
492 MFillH fill0(&mh1);
493 //fill0.SetDrawOption("projx only");
494
495 MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
496 MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
497 if (bins2 && bins3)
498 {
499 bins2->SetName("BinningThetaEY");
500 bins3->SetName("BinningThetaEX");
501 }
502 tlist1.AddToList(&readmc);
503 tlist1.AddToList(&weight);
504
505 temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
506 MH3 mh3mc(temp1);
507
508 MFEventSelector2 sel1(mh3mc);
509 sel1.SetHistIsProbability();
510
511 fill0.SetFilter(&sel1);
512
513 if (!fRawMc)
514 tlist1.AddToList(&sel1);
515 tlist1.AddToList(&fill0);
516
517 MEvtLoop loop1(fName);
518 loop1.SetParList(&plist);
519 loop1.SetLogStream(fLog);
520 loop1.SetDisplay(fDisplay);
521
522 if (!SetupEnv(loop1))
523 return kFALSE;
524
525 if (!loop1.Eventloop(fMaxEvents))
526 {
527 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
528 return kFALSE;
529 }
530
531 if (!loop1.GetDisplay())
532 {
533 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
534 return kFALSE;
535 }
536
537 if (bins2 && bins3)
538 {
539 bins2->SetName("BinningEnergyEst");
540 bins3->SetName("BinningTheta");
541 }
542
543 return kTRUE;
544}
545
546// --------------------------------------------------------------------------
547//
548// Calculate the final spectrum from:
549// - collection area
550// - excess
551// - correction coefficients
552// - ontime
553// and display it
554//
555TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
556{
557 TH1D collarea(area.GetHEnergy());
558 TH1D spectrum(excess);
559 TH1D weights;
560 hest.GetWeights(weights);
561
562 cout << "Effective On time: " << ontime << "s" << endl;
563
564 spectrum.SetDirectory(NULL);
565 spectrum.SetBit(kCanDelete);
566 spectrum.Scale(1./ontime);
567 spectrum.Divide(&collarea);
568 spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
569 spectrum.SetYTitle("N/sm^{2}");
570
571 TCanvas &c1 = fDisplay->AddTab("Spectrum");
572 c1.Divide(2,2);
573 c1.cd(1);
574 gPad->SetBorderMode(0);
575 gPad->SetLogx();
576 gPad->SetLogy();
577 gPad->SetGridx();
578 gPad->SetGridy();
579 collarea.DrawCopy();
580
581 c1.cd(2);
582 gPad->SetBorderMode(0);
583 gPad->SetLogx();
584 gPad->SetLogy();
585 gPad->SetGridx();
586 gPad->SetGridy();
587 spectrum.DrawCopy();
588
589 c1.cd(3);
590 gPad->SetBorderMode(0);
591 gPad->SetLogx();
592 gPad->SetLogy();
593 gPad->SetGridx();
594 gPad->SetGridy();
595 weights.DrawCopy();
596
597 //spectrum.Multiply(&weights);
598 spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
599 spectrum.SetBit(TH1::kNoStats);
600
601 for (int i=0; i<excess.GetNbinsX(); i++)
602 {
603 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
604 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1));
605
606 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
607 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
608 }
609
610 c1.cd(4);
611 gPad->SetBorderMode(0);
612 gPad->SetLogx();
613 gPad->SetLogy();
614 gPad->SetGridx();
615 gPad->SetGridy();
616 spectrum.SetMinimum(1e-12);
617 spectrum.SetXTitle("E [GeV]");
618 spectrum.SetYTitle("N/TeVsm^{2}");
619 spectrum.DrawCopy();
620
621 TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
622 f.SetParameter(0, -2.87);
623 f.SetParameter(1, 1.9e-6);
624 f.SetLineColor(kGreen);
625 spectrum.Fit(&f, "NIM", "", 100, 4000);
626 f.DrawCopy("same");
627
628 const Double_t p0 = f.GetParameter(0);
629 const Double_t p1 = f.GetParameter(1);
630
631 const Double_t e0 = f.GetParError(0);
632 const Double_t e1 = f.GetParError(1);
633
634 const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
635 const Double_t exp = TMath::Power(10, np);
636
637 TString str;
638 str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
639 str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
640 str += "\\frac{ph}{TeVm^{2}s}";
641
642 TLatex tex;
643 tex.SetTextSize(0.045);
644 tex.SetBit(TLatex::kTextNDC);
645 tex.DrawLatex(0.45, 0.935, str);
646
647 str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
648 tex.DrawLatex(0.70, 0.83, str);
649
650 TArrayD res(2);
651 res[0] = f.GetParameter(0);
652 res[1] = f.GetParameter(1);
653
654/*
655 // Plot other spectra from Whipple
656 f.SetParameter(0, -2.45);
657 f.SetParameter(1, 3.3e-7);
658 f.SetRange(300, 8000);
659 f.SetLineColor(kBlack);
660 f.SetLineStyle(kDashed);
661 f.DrawCopy("same");
662
663 // Plot other spectra from Cangaroo
664 f.SetParameter(0, -2.53);
665 f.SetParameter(1, 2.0e-7);
666 f.SetRange(7000, 50000);
667 f.SetLineColor(kBlack);
668 f.SetLineStyle(kDashed);
669 f.DrawCopy("same");
670
671 // Plot other spectra from Robert
672 f.SetParameter(0, -2.59);
673 f.SetParameter(1, 2.58e-7);
674 f.SetRange(150, 1500);
675 f.SetLineColor(kBlack);
676 f.SetLineStyle(kDashed);
677 f.DrawCopy("same");
678
679 // Plot other spectra from HEGRA
680 f.SetParameter(0, -2.61);
681 f.SetParameter(1, 2.7e-7);
682 f.SetRange(1000, 20000);
683 f.SetLineColor(kBlack);
684 f.SetLineStyle(kDashed);
685 f.DrawCopy("same");
686 */
687 return res;
688}
689
690// --------------------------------------------------------------------------
691//
692// Scale some image parameter plots using the scale factor and plot them
693// together with the corresponding MC histograms.
694// Called from DisplaySize
695//
696Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
697{
698 TString same(name);
699 same += "Same";
700
701 TH1 *h1 = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
702 TH1 *h2 = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
703 if (!h1 || !h2)
704 return kFALSE;
705
706 TObject *obj = plist.FindObject(plot);
707 if (!obj)
708 {
709 *fLog << warn << plot << " not in parameter list... skipping." << endl;
710 return kFALSE;
711 }
712
713 TH1 *h3 = (TH1*)obj->FindObject(name);
714 if (!h3)
715 {
716 *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
717 return kFALSE;
718 }
719
720
721 const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
722 const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
723
724 gPad->SetBorderMode(0);
725 h2->SetLineColor(kBlack);
726 h3->SetLineColor(kBlue);
727 h2->Add(h1, -ascale);
728
729 //h2->Scale(1./ontime); //h2->Integral());
730 h3->Scale(scale); //h3->Integral());
731
732 h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
733
734 h2 = h2->DrawCopy();
735 h3 = h3->DrawCopy("same");
736
737 // Don't do this on the original object!
738 h2->SetStats(kFALSE);
739 h3->SetStats(kFALSE);
740
741 return kTRUE;
742}
743
744// --------------------------------------------------------------------------
745//
746// Take a lot of histograms and plot them together in one plot.
747// Calls PlotSame
748//
749Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
750{
751 *fLog << inf << "Reading from file: " << fPathIn << endl;
752
753 TFile file(fPathIn, "READ");
754 if (!file.IsOpen())
755 {
756 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
757 return kFALSE;
758 }
759
760 file.cd();
761 MStatusArray arr;
762 if (arr.Read()<=0)
763 {
764 *fLog << "MStatusDisplay not found in file... abort." << endl;
765 return kFALSE;
766 }
767
768 TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
769 if (!excess)
770 return kFALSE;
771
772 // ------------------- Plot excess versus size -------------------
773
774 TCanvas &c = fDisplay->AddTab("Excess");
775 c.Divide(3,2);
776 c.cd(1);
777 gPad->SetBorderMode(0);
778 gPad->SetLogx();
779 gPad->SetLogy();
780 gPad->SetGridx();
781 gPad->SetGridy();
782
783 excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
784 excess = excess->DrawCopy("E2");
785 // Don't do this on the original object!
786 excess->SetStats(kFALSE);
787 excess->SetMarkerStyle(kFullDotMedium);
788 excess->SetFillColor(kBlack);
789 excess->SetFillStyle(0);
790 excess->SetName("Excess ");
791 excess->SetDirectory(0);
792
793 TObject *o=0;
794 if ((o=plist.FindObject("ExcessMC")))
795 {
796 TH1 *histsel = (TH1F*)o->FindObject("");
797 if (histsel)
798 {
799 if (scale<0)
800 scale = excess->Integral()/histsel->Integral();
801
802 histsel->Scale(scale);
803 histsel->SetLineColor(kBlue);
804 histsel->SetBit(kCanDelete);
805 histsel = histsel->DrawCopy("E1 same");
806 // Don't do this on the original object!
807 histsel->SetStats(kFALSE);
808
809 fLog->Separator("Kolmogorov Test");
810 histsel->KolmogorovTest(excess, "DX");
811 fLog->Separator("Chi^2 Test");
812 const Double_t p = histsel->Chi2Test(excess, "P");
813
814 TLatex tex;
815 tex.SetBit(TLatex::kTextNDC);
816 tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
817 }
818 }
819
820 // -------------- Comparison of Image Parameters --------------
821 c.cd(2);
822 PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale);
823
824 c.cd(3);
825 PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
826
827 c.cd(4);
828 PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale);
829
830 c.cd(5);
831 PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale);
832
833 c.cd(6);
834 PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale);
835
836 return kTRUE;
837}
838
839Bool_t MJSpectrum::Process(const MDataSet &set)
840{
841 if (!set.IsValid())
842 {
843 *fLog << err << "ERROR - DataSet invalid!" << endl;
844 return kFALSE;
845 }
846
847 CheckEnv();
848
849 // --------------------------------------------------------------------------------
850
851 *fLog << inf;
852 fLog->Separator(GetDescriptor());
853 *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
854 *fLog << endl;
855
856 MBinning bins1("BinningAlpha");
857 MBinning bins2("BinningEnergyEst");
858 MBinning bins3("BinningTheta");
859 MBinning bins4("BinningFalseSource");
860 MBinning bins5("BinningWidth");
861 MBinning bins6("BinningLength");
862 MBinning bins7("BinningDist");
863 MBinning bins8("BinningMaxDist");
864 MBinning bins9("BinningM3Long");
865 MBinning bins0("BinningConc1");
866
867 MAlphaFitter fit;
868
869 MParList plist;
870 plist.AddToList(&bins1);
871 plist.AddToList(&bins2);
872 plist.AddToList(&bins3);
873 plist.AddToList(&bins4);
874 plist.AddToList(&bins5);
875 plist.AddToList(&bins6);
876 plist.AddToList(&bins7);
877 plist.AddToList(&bins8);
878 plist.AddToList(&bins9);
879 plist.AddToList(&bins0);
880 plist.AddToList(&fit);
881
882 TH1D temp1, size;
883 const Float_t ontime = ReadInput(plist, temp1, size);
884 if (ontime<0)
885 {
886 *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
887 return kFALSE;
888 }
889
890 MMcSpectrumWeight weight;
891 if (!InitWeighting(set, weight))
892 return kFALSE;
893
894 PrintSetup(fit);
895 bins3.SetEdges(temp1, 'x');
896
897 TH1D temp2(temp1);
898 if (!ReadOrigMCDistribution(set, temp2, weight))
899 return kFALSE;
900
901 if (!GetThetaDistribution(temp1, temp2))
902 return kFALSE;
903
904 if (!fNoThetaWeights)
905 weight.SetZdWeights(&temp1);
906
907 TH1D excess;
908 if (!Refill(plist, excess))
909 return kFALSE;
910
911 TH2D hist;
912 MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
913 if (fSimpleMode)
914 {
915 hist.UseCurrentStyle();
916 MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
917 if (!ReadOrigMCDistribution(set, hist, weight))
918 return kFALSE;
919
920 if (!fRawMc)
921 {
922 for (int y=0; y<hist.GetNbinsY(); y++)
923 for (int x=0; x<hist.GetNbinsX(); x++)
924 hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
925 //hist.SetEntries(hist.Integral());
926 }
927 }
928 else
929 {
930 weight.SetNameMcEvt("MMcEvtBasic");
931 if (!IntermediateLoop(plist, mh1, temp1, set, weight))
932 return kFALSE;
933 weight.SetNameMcEvt();
934 }
935
936 DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
937
938 // ------------------------- Final loop --------------------------
939
940 *fLog << endl;
941 fLog->Separator("Calculate efficiencies");
942 *fLog << endl;
943
944 MTaskList tlist2;
945 plist.Replace(&tlist2);
946
947 MReadMarsFile read("Events");
948 read.DisableAutoScheme();
949 set.AddFilesOn(read);
950
951 // Selector to get correct (final) theta-distribution
952 temp1.SetXTitle("MPointingPos.fZd");
953
954 MH3 mh3(temp1);
955
956 MFEventSelector2 sel2(mh3);
957 sel2.SetHistIsProbability();
958
959 MContinue contsel(&sel2);
960 contsel.SetInverted();
961
962 // Get correct source position
963 MSrcPosCalc calc;
964
965 // Calculate corresponding Hillas parameters
966 MHillasCalc hcalc1;
967 MHillasCalc hcalc2("MHillasCalcAnti");
968 hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
969 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
970 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
971 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
972
973 // Fill collection area and energy estimator (unfolding)
974 // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
975 MHCollectionArea area;
976 area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
977 MHEnergyEst hest;
978
979 MFillH fill3(&area, "", "FillCollectionArea");
980 MFillH fill4(&hest, "", "FillEnergyEst");
981 fill3.SetWeight();
982 fill4.SetWeight();
983
984 MH3 hsize("MHillas.fSize");
985 hsize.SetName("ExcessMC");
986 hsize.Sumw2();
987
988 MBinning bins(size, "BinningExcessMC");
989 plist.AddToList(&hsize);
990 plist.AddToList(&bins);
991
992 MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
993 MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
994 MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
995 MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
996 MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
997 MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
998 MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
999 MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC");
1000 fill1a.SetNameTab("PreCut");
1001 fill2a.SetNameTab("PostCut");
1002 fill3a.SetNameTab("VsSize");
1003 fill4a.SetNameTab("HilExt");
1004 fill5a.SetNameTab("HilSrc");
1005 fill6a.SetNameTab("ImgPar");
1006 fill7a.SetNameTab("NewPar");
1007 fill8a.SetBit(MFillH::kDoNotDisplay);
1008 fill1a.SetWeight();
1009 fill2a.SetWeight();
1010 fill3a.SetWeight();
1011 fill4a.SetWeight();
1012 fill5a.SetWeight();
1013 fill6a.SetWeight();
1014 fill7a.SetWeight();
1015 fill8a.SetWeight();
1016
1017 MEnergyEstimate est;
1018 MTaskEnv taskenv1("EstimateEnergy");
1019 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
1020
1021 tlist2.AddToList(&read);
1022 if (!fRawMc && fNoThetaWeights)
1023 tlist2.AddToList(&contsel);
1024 tlist2.AddToList(&calc);
1025 tlist2.AddToList(&hcalc1);
1026 tlist2.AddToList(&hcalc2);
1027 tlist2.AddToList(&weight);
1028 tlist2.AddToList(&fill1a);
1029 tlist2.AddToList(fCut0);
1030 tlist2.AddToList(fCut1);
1031 tlist2.AddToList(fCut2);
1032 tlist2.AddToList(fCut3);
1033 tlist2.AddToList(&taskenv1);
1034 tlist2.AddToList(&fill3);
1035 tlist2.AddToList(&fill4);
1036 tlist2.AddToList(&fill2a);
1037 tlist2.AddToList(&fill3a);
1038 tlist2.AddToList(&fill4a);
1039 tlist2.AddToList(&fill5a);
1040 tlist2.AddToList(&fill6a);
1041 tlist2.AddToList(&fill7a);
1042 tlist2.AddToList(&fill8a);
1043 //tlist2.AddToList(&fill9a);
1044
1045 MEvtLoop loop2(fName);
1046 loop2.SetParList(&plist);
1047 loop2.SetDisplay(fDisplay);
1048 loop2.SetLogStream(fLog);
1049
1050 if (!SetupEnv(loop2))
1051 return kFALSE;
1052
1053 if (!loop2.Eventloop(fMaxEvents))
1054 {
1055 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
1056 return kFALSE;
1057 }
1058
1059 if (!loop2.GetDisplay())
1060 {
1061 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
1062 return kFALSE;
1063 }
1064
1065 gLog.Separator("Energy Estimator");
1066 if (plist.FindObject("EstimateEnergy"))
1067 plist.FindObject("EstimateEnergy")->Print();
1068
1069 gLog.Separator("Spectrum");
1070
1071 // -------------------------- Spectrum ----------------------------
1072
1073 // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
1074 TArrayD res(DisplaySpectrum(area, excess, hest, ontime));
1075
1076 // Spectrum fitted (convert res[1] from TeV to GeV)
1077 TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
1078
1079 // Number of events this spectrum would produce per s and m^2
1080 Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
1081
1082 // scale with effective collection area to get the event rate (N/s)
1083 // scale with the effective observation time to absolute observed events
1084 n *= area.GetCollectionAreaAbs()*ontime; // N
1085
1086 // Now calculate the scale factor from the number of events
1087 // produced and the number of events which should have been
1088 // observed with our telescope in the time ontime
1089 const Double_t scale = n/area.GetEntries();
1090
1091 // Print normalization constant
1092 cout << "MC normalization factor: " << scale << endl;
1093
1094 // Overlay normalized plots
1095 DisplaySize(plist, scale);
1096
1097 // check if output should be written
1098 if (!fPathOut.IsNull())
1099 fDisplay->SaveAsRoot(fPathOut);
1100
1101 return kTRUE;
1102}
Note: See TracBrowser for help on using the repository browser.