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

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