source: tags/Mars-V0.9.2/mjobs/MJSpectrum.cc

Last change on this file was 7071, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 23.4 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 <TCanvas.h>
41#include <TObjArray.h>
42
43// Environment
44#include "MLog.h"
45#include "MLogManip.h"
46
47#include "MStatusArray.h"
48#include "MStatusDisplay.h"
49
50// Container
51#include "MH3.h"
52#include "MBinning.h"
53#include "MDataSet.h"
54
55// Spectrum
56#include "../mhflux/MAlphaFitter.h"
57#include "../mhflux/MHAlpha.h"
58#include "../mhflux/MHCollectionArea.h"
59#include "../mhflux/MHEnergyEst.h"
60
61// Eventloop
62#include "MEvtLoop.h"
63#include "MTaskList.h"
64#include "MParList.h"
65
66// Tasks/Filter
67#include "MReadMarsFile.h"
68#include "MReadMarsFile.h"
69#include "MFEventSelector2.h"
70#include "MFDataMember.h"
71#include "MEnergyEstimate.h"
72#include "MTaskEnv.h"
73#include "MFillH.h"
74#include "MHillasCalc.h"
75#include "MSrcPosCalc.h"
76#include "MContinue.h"
77
78ClassImp(MJSpectrum);
79
80using namespace std;
81
82MJSpectrum::MJSpectrum(const char *name, const char *title)
83 : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
84 fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE)
85{
86 fName = name ? name : "MJSpectrum";
87 fTitle = title ? title : "Standard program to calculate spectrum";
88}
89
90MJSpectrum::~MJSpectrum()
91{
92 if (fCut0)
93 delete fCut0;
94 if (fCut1)
95 delete fCut1;
96 if (fCut2)
97 delete fCut2;
98 if (fCut3)
99 delete fCut3;
100 if (fEstimateEnergy)
101 delete fEstimateEnergy;
102}
103
104// --------------------------------------------------------------------------
105//
106// Setup a task estimating the energy. The given task is cloned.
107//
108void MJSpectrum::SetEnergyEstimator(const MTask *task)
109{
110 if (fEstimateEnergy)
111 delete fEstimateEnergy;
112 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
113}
114
115Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const
116{
117 if (task)
118 {
119 delete task;
120 task = 0;
121 }
122
123 task = (MTask*)gFile->Get(name);
124 if (!task)
125 {
126 *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
127 return kFALSE;
128 }
129 if (!task->InheritsFrom(MTask::Class()))
130 {
131 *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
132 delete task;
133 return kFALSE;
134 }
135
136 task->SetName(name);
137
138 return kTRUE;
139}
140
141void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
142{
143 fLog->Separator("Alpha Fitter");
144 *fLog << all;
145 fit.Print();
146
147 fLog->Separator("Used Cuts");
148 fCut0->Print();
149 fCut1->Print();
150 fCut2->Print();
151 fCut3->Print();
152
153 //gLog.Separator("Energy Estimator");
154 //fEstimateEnergy->Print();
155}
156
157Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
158{
159 *fLog << inf << "Reading from file: " << fPathIn << endl;
160
161 TFile file(fPathIn, "READ");
162 if (!file.IsOpen())
163 {
164 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
165 return -1;
166 }
167
168 MStatusArray arr;
169 if (arr.Read()<=0)
170 {
171 *fLog << "MStatusDisplay not found in file... abort." << endl;
172 return -1;
173 }
174
175 TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
176 TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
177 if (!vstime || !size)
178 return -1;
179
180 vstime->Copy(h1);
181 size->Copy(h2);
182 h1.SetDirectory(0);
183 h2.SetDirectory(0);
184
185 if (fDisplay)
186 arr.DisplayIn(*fDisplay, "Hist");
187
188 if (!ReadTask(fCut0, "Cut0"))
189 return -1;
190 if (!ReadTask(fCut1, "Cut1"))
191 return -1;
192 if (!ReadTask(fCut2, "Cut2"))
193 return -1;
194 if (!ReadTask(fCut3, "Cut3"))
195 return -1;
196
197 TObjArray arrread;
198
199 TIter Next(plist);
200 TObject *o=0;
201 while ((o=Next()))
202 if (o->InheritsFrom(MBinning::Class()))
203 arrread.Add(o);
204
205 arrread.Add(plist.FindObject("MAlphaFitter"));
206
207 if (!ReadContainer(arrread))
208 return -1;
209
210 return vstime->Integral();
211}
212
213Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const
214{
215 // Some debug output
216 fLog->Separator("Compiling original MC distribution");
217
218 *fLog << inf << "Please stand by, this may take a while..." << flush;
219
220 if (fDisplay)
221 fDisplay->SetStatusLine1("Compiling MC distribution...");
222
223 // Create chain
224 TChain chain("OriginalMC");
225 set.AddFilesOn(chain);
226
227 // Prepare histogram
228 h.Reset();
229
230 // Fill histogram from chain
231 h.SetDirectory(gROOT);
232 if (h.InheritsFrom(TH2::Class()))
233 {
234 h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
235 h.SetXTitle("\\Theta [\\circ]");
236 h.SetYTitle("E [GeV]");
237 h.SetZTitle("Counts");
238 chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff");
239 }
240 else
241 {
242 h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
243 h.SetXTitle("\\Theta [\\circ]");
244 h.SetYTitle("Counts");
245 chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff");
246 }
247 h.SetDirectory(0);
248
249 *fLog << "done." << endl;
250 if (fDisplay)
251 fDisplay->SetStatusLine2("done.");
252
253 if (h.GetEntries()>0)
254 return kTRUE;
255
256 *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
257
258 return h.GetEntries()>0;
259}
260
261Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
262{
263 // Display some stuff
264 if (fDisplay)
265 {
266 TCanvas &c = fDisplay->AddTab("ZdDist");
267 c.Divide(2,2);
268
269 // On-Time vs. Theta
270 c.cd(1);
271 gPad->SetBorderMode(0);
272 temp1.DrawCopy();
273
274 // Number of MC events (produced) vs Theta
275 c.cd(2);
276 gPad->SetBorderMode(0);
277 temp2.SetName("NVsTheta");
278 temp2.DrawCopy();
279
280 c.cd(4);
281 gPad->SetBorderMode(0);
282
283 c.cd(3);
284 gPad->SetBorderMode(0);
285 }
286
287 // Calculate the Probability
288 temp1.Divide(&temp2);
289 temp1.Scale(1./temp1.GetMaximum());
290
291 // Some cosmetics: Name, Axis, etc.
292 temp1.SetName("ProbVsTheta");
293 temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
294 temp1.SetYTitle("Probability");
295 if (fDisplay)
296 temp1.DrawCopy();
297
298 return kTRUE;
299}
300
301void MJSpectrum::DisplayResult(const TH2D &h2) const
302{
303 if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
304 return;
305
306 TH1D &proj = *h2.ProjectionX();
307 proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
308 proj.SetXTitle("\\Theta [\\circ]");
309 proj.SetYTitle("Counts");
310 proj.SetLineColor(kBlue);
311 proj.SetDirectory(0);
312 proj.SetBit(kCanDelete);
313
314 TVirtualPad *pad = gPad;
315
316 pad->cd(4);
317 proj.DrawCopy();
318
319 pad->cd(1);
320 TH1D *theta = (TH1D*)gPad->FindObject("Theta");
321 if (theta)
322 {
323 proj.Scale(theta->GetMaximum()/proj.GetMaximum());
324 theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
325 }
326 proj.Draw("same");
327}
328
329Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
330{
331 // Try to find the class used to determin the signal!
332 TString cls("MHAlpha");
333 if (fDisplay)
334 {
335 TCanvas *c = fDisplay->GetCanvas("Hist");
336 if (c)
337 {
338 TIter Next(c->GetListOfPrimitives());
339 TObject *obj=0;
340 while ((obj=Next()))
341 if (obj->InheritsFrom(MHAlpha::Class()))
342 break;
343 if (obj)
344 cls = obj->ClassName();
345 }
346 }
347
348 cout << "FOUND: "<< cls << endl;
349
350 // Now fill the histogram
351 *fLog << endl;
352 fLog->Separator("Refill Excess");
353 *fLog << endl;
354
355 MTaskList tlist;
356 plist.AddToList(&tlist);
357
358 MReadTree read("Events");
359 read.DisableAutoScheme();
360 read.AddFile(fPathIn);
361
362 MEnergyEstimate est;
363 MTaskEnv taskenv1("EstimateEnergy");
364 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
365
366 MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
367 MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
368
369 MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
370 MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
371
372 fill1.SetFilter(&f0);
373 fill2.SetFilter(&f1);
374
375 tlist.AddToList(&read);
376 tlist.AddToList(&taskenv1);
377 tlist.AddToList(&f0);
378 tlist.AddToList(&f1);
379 tlist.AddToList(&fill1);
380 tlist.AddToList(&fill2);
381
382 MEvtLoop loop(fName);
383 loop.SetParList(&plist);
384 loop.SetDisplay(fDisplay);
385 loop.SetLogStream(fLog);
386
387 if (!SetupEnv(loop))
388 return kFALSE;
389
390 if (!loop.Eventloop())
391 {
392 *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
393 return kFALSE;
394 }
395
396 if (!loop.GetDisplay())
397 {
398 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
399 return kFALSE;
400 }
401
402 const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE");
403 if (!halpha)
404 {
405 *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
406 return kFALSE;
407 }
408
409 halpha->GetHEnergy().Copy(h2);
410 h2.SetDirectory(0);
411
412 return kTRUE;
413}
414
415Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
416{
417 MTaskList tlist1;
418 plist.Replace(&tlist1);
419
420 MReadMarsFile readmc("OriginalMC");
421 //readmc.DisableAutoScheme();
422 set.AddFilesOn(readmc);
423 readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
424 readmc.EnableBranch("MMcEvtBasic.fEnergy");
425
426 mh1.SetLogy();
427 mh1.SetLogz();
428 mh1.SetName("ThetaE");
429
430 MFillH fill0(&mh1);
431 //fill0.SetDrawOption("projx only");
432
433 MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
434 MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
435 if (bins2 && bins3)
436 {
437 bins2->SetName("BinningThetaEY");
438 bins3->SetName("BinningThetaEX");
439 }
440 tlist1.AddToList(&readmc);
441
442 temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
443 MH3 mh3mc(temp1);
444
445 MFEventSelector2 sel1(mh3mc);
446 sel1.SetHistIsProbability();
447
448 fill0.SetFilter(&sel1);
449
450 if (!fRawMc)
451 tlist1.AddToList(&sel1);
452 tlist1.AddToList(&fill0);
453
454 MEvtLoop loop1(fName);
455 loop1.SetParList(&plist);
456 loop1.SetLogStream(fLog);
457 loop1.SetDisplay(fDisplay);
458
459 if (!SetupEnv(loop1))
460 return kFALSE;
461
462 if (!loop1.Eventloop(fMaxEvents))
463 {
464 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
465 return kFALSE;
466 }
467
468 if (!loop1.GetDisplay())
469 {
470 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
471 return kFALSE;
472 }
473
474 if (bins2 && bins3)
475 {
476 bins2->SetName("BinningEnergyEst");
477 bins3->SetName("BinningTheta");
478 }
479 return kTRUE;
480}
481
482void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
483{
484 TH1D collarea(area.GetHEnergy());
485 TH1D spectrum(excess);
486 TH1D weights;
487 hest.GetWeights(weights);
488
489 cout << "Effective On time: " << ontime << "s" << endl;
490
491 spectrum.SetDirectory(NULL);
492 spectrum.SetBit(kCanDelete);
493 spectrum.Scale(1./ontime);
494 spectrum.Divide(&collarea);
495 spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
496 spectrum.SetYTitle("N/sm^{2}");
497
498 TCanvas &c1 = fDisplay->AddTab("Spectrum");
499 c1.Divide(2,2);
500 c1.cd(1);
501 gPad->SetBorderMode(0);
502 gPad->SetLogx();
503 gPad->SetLogy();
504 gPad->SetGridx();
505 gPad->SetGridy();
506 collarea.DrawCopy();
507
508 c1.cd(2);
509 gPad->SetBorderMode(0);
510 gPad->SetLogx();
511 gPad->SetLogy();
512 gPad->SetGridx();
513 gPad->SetGridy();
514 spectrum.DrawCopy();
515
516 c1.cd(3);
517 gPad->SetBorderMode(0);
518 gPad->SetLogx();
519 gPad->SetLogy();
520 gPad->SetGridx();
521 gPad->SetGridy();
522 weights.DrawCopy();
523
524 //spectrum.Divide(&weights);
525 spectrum.Multiply(&weights);
526 spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
527
528 for (int i=0; i<excess.GetNbinsX(); i++)
529 {
530 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
531 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
532 }
533
534 c1.cd(4);
535 gPad->SetBorderMode(0);
536 gPad->SetLogx();
537 gPad->SetLogy();
538 gPad->SetGridx();
539 gPad->SetGridy();
540 spectrum.SetXTitle("E [GeV]");
541 spectrum.SetYTitle("N/TeVsm^{2}");
542 spectrum.DrawCopy();
543
544 TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
545 f.SetParameter(0, -2.87);
546 f.SetParameter(1, 1.9e-6);
547 f.SetLineColor(kGreen);
548 spectrum.Fit(&f, "NI", "", 55, 2e4);
549 f.DrawCopy("same");
550
551 /*
552 TString str;
553 str += "(1.68#pm0.15)10^{-7}";
554 str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
555 str += "\\frac{ph}{TeVm^{2}s}";
556
557 TLatex tex;
558 tex.DrawLatex(2e2, 7e-5, str);
559 */
560}
561
562Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const
563{
564 cout << name << endl;
565 TString same(name);
566 same += "Same";
567
568 TH1 *h1 = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
569 TH1 *h2 = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
570 if (!h1 || !h2)
571 return kFALSE;
572
573 TObject *obj = plist.FindObject(plot);
574 if (!obj)
575 {
576 *fLog << warn << plot << " not in parameter list... skipping." << endl;
577 return kFALSE;
578 }
579
580 TH1 *h3 = (TH1*)obj->FindObject(name);
581 if (!h3)
582 {
583 *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
584 return kFALSE;
585 }
586
587
588 const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
589 const Double_t scale = fit ? fit->GetScaleFactor() : 1;
590
591 gPad->SetBorderMode(0);
592 h2->SetLineColor(kBlack);
593 h3->SetLineColor(kBlue);
594 h2->Add(h1, -scale);
595
596 h2->Scale(1./h2->Integral());
597 h3->Scale(1./h3->Integral());
598
599 h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
600
601 h2 = h2->DrawCopy();
602 h3 = h3->DrawCopy("same");
603
604 // Don't do this on the original object!
605 h2->SetStats(kFALSE);
606 h3->SetStats(kFALSE);
607
608 return kTRUE;
609}
610
611Bool_t MJSpectrum::DisplaySize(MParList &plist) const
612{
613 *fLog << inf << "Reading from file: " << fPathIn << endl;
614
615 TFile file(fPathIn, "READ");
616 if (!file.IsOpen())
617 {
618 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
619 return kFALSE;
620 }
621
622 file.cd();
623 MStatusArray arr;
624 if (arr.Read()<=0)
625 {
626 *fLog << "MStatusDisplay not found in file... abort." << endl;
627 return kFALSE;
628 }
629
630 TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
631 if (!excess)
632 return kFALSE;
633
634 // ------------------- Plot excess versus size -------------------
635
636 TCanvas &c = fDisplay->AddTab("Excess");
637 c.Divide(3,2);
638 c.cd(1);
639 gPad->SetBorderMode(0);
640 gPad->SetLogx();
641 gPad->SetLogy();
642 gPad->SetGridx();
643 gPad->SetGridy();
644
645 excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
646 excess->Scale(1./excess->Integral());
647 excess = excess->DrawCopy();
648 // Don't do this on the original object!
649 excess->SetStats(kFALSE);
650
651 TObject *o=0;
652 if ((o=plist.FindObject("ExcessSize")))
653 {
654 TH1 *histsel = (TH1F*)o->FindObject("");
655 if (histsel)
656 {
657 histsel->Scale(1./histsel->Integral());
658 histsel->SetLineColor(kBlue);
659 histsel->SetBit(kCanDelete);
660 histsel = histsel->DrawCopy("same");
661 // Don't do this on the original object!
662 histsel->SetStats(kFALSE);
663 }
664 }
665
666 // -------------- Comparison of Image Parameters --------------
667 c.cd(2);
668 PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost");
669
670 c.cd(3);
671 PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
672
673 c.cd(4);
674 PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost");
675
676 c.cd(5);
677 PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost");
678
679 c.cd(6);
680 PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost");
681
682 return kTRUE;
683}
684
685Bool_t MJSpectrum::Process(const MDataSet &set)
686{
687 if (!set.IsValid())
688 {
689 *fLog << err << "ERROR - DataSet invalid!" << endl;
690 return kFALSE;
691 }
692
693 CheckEnv();
694
695 // --------------------------------------------------------------------------------
696
697 *fLog << inf;
698 fLog->Separator(GetDescriptor());
699 *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
700 *fLog << endl;
701
702 MBinning bins1("BinningAlpha");
703 MBinning bins2("BinningEnergyEst");
704 MBinning bins3("BinningTheta");
705 MBinning bins4("BinningFalseSource");
706 MBinning bins5("BinningWidth");
707 MBinning bins6("BinningLength");
708 MBinning bins7("BinningDist");
709 MBinning bins8("BinningMaxDist");
710 MBinning bins9("BinningM3Long");
711 MBinning bins0("BinningConc1");
712
713 MAlphaFitter fit;
714
715 MParList plist;
716 plist.AddToList(&bins1);
717 plist.AddToList(&bins2);
718 plist.AddToList(&bins3);
719 plist.AddToList(&bins4);
720 plist.AddToList(&bins5);
721 plist.AddToList(&bins6);
722 plist.AddToList(&bins7);
723 plist.AddToList(&bins8);
724 plist.AddToList(&bins9);
725 plist.AddToList(&bins0);
726 plist.AddToList(&fit);
727
728 TH1D temp1, size;
729 const Float_t ontime = ReadInput(plist, temp1, size);
730 if (ontime<0)
731 {
732 *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
733 return kFALSE;
734 }
735
736 PrintSetup(fit);
737 bins3.SetEdges(temp1, 'x');
738
739 TH1D temp2(temp1);
740 if (!ReadOrigMCDistribution(set, temp2))
741 return kFALSE;
742
743 if (!GetThetaDistribution(temp1, temp2))
744 return kFALSE;
745
746 TH1D excess;
747 if (!Refill(plist, excess))
748 return kFALSE;
749
750 TH2D hist;
751 MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
752 if (fSimpleMode)
753 {
754 hist.UseCurrentStyle();
755 MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
756 if (!ReadOrigMCDistribution(set, hist))
757 return kFALSE;
758
759 if (!fRawMc)
760 {
761 for (int y=0; y<hist.GetNbinsY(); y++)
762 for (int x=0; x<hist.GetNbinsX(); x++)
763 hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
764 }
765 }
766 else
767 if (!IntermediateLoop(plist, mh1, temp1, set))
768 return kFALSE;
769
770 DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
771
772 // ------------------------- Final loop --------------------------
773
774 *fLog << endl;
775 fLog->Separator("Calculate efficiencies");
776 *fLog << endl;
777
778 MTaskList tlist2;
779 plist.Replace(&tlist2);
780
781 MReadMarsFile read("Events");
782 read.DisableAutoScheme();
783 set.AddFilesOn(read);
784
785 // Selector to get correct (final) theta-distribution
786 temp1.SetXTitle("MPointingPos.fZd");
787
788 MH3 mh3(temp1);
789
790 MFEventSelector2 sel2(mh3);
791 sel2.SetHistIsProbability();
792
793 MContinue contsel(&sel2);
794 contsel.SetInverted();
795
796 // Get correct source position
797 MSrcPosCalc calc;
798
799 // Calculate corresponding Hillas parameters
800 MHillasCalc hcalc1;
801 MHillasCalc hcalc2("MHillasCalcAnti");
802 hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
803 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
804 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
805 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
806
807 // Fill collection area and energy estimator (unfolding)
808 // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
809 MHCollectionArea area;
810 area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
811 MHEnergyEst hest;
812
813 MFillH fill3(&area, "", "FillCollectionArea");
814 MFillH fill4(&hest, "", "FillEnergyEst");
815
816 MH3 hsize("MHillas.fSize");
817 //MH3 henergy("MEnergyEst.fVal");
818 hsize.SetName("ExcessSize");
819 //henergy.SetName("EnergyEst");
820 MBinning bins(size, "BinningExcessSize");
821 plist.AddToList(&hsize);
822 //plist.AddToList(&henergy);
823 plist.AddToList(&bins);
824
825 MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
826 MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
827 MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
828 MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
829 MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
830 MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
831 MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
832 MFillH fill8a("ExcessSize [MH3]", "", "FillExcessSize");
833 //MFillH fill9a("EnergyEst [MH3]", "", "FillExcessEEst");
834 fill1a.SetNameTab("PreCut");
835 fill2a.SetNameTab("PostCut");
836 fill3a.SetNameTab("VsSize");
837 fill4a.SetNameTab("HilExt");
838 fill5a.SetNameTab("HilSrc");
839 fill6a.SetNameTab("ImgPar");
840 fill7a.SetNameTab("NewPar");
841 fill8a.SetBit(MFillH::kDoNotDisplay);
842 //fill9a.SetBit(MFillH::kDoNotDisplay);
843
844 MEnergyEstimate est;
845 MTaskEnv taskenv1("EstimateEnergy");
846 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
847
848 tlist2.AddToList(&read);
849 if (!fRawMc)
850 tlist2.AddToList(&contsel);
851 tlist2.AddToList(&calc);
852 tlist2.AddToList(&hcalc1);
853 tlist2.AddToList(&hcalc2);
854 tlist2.AddToList(&fill1a);
855 tlist2.AddToList(fCut0);
856 tlist2.AddToList(fCut1);
857 tlist2.AddToList(fCut2);
858 tlist2.AddToList(fCut3);
859 tlist2.AddToList(&taskenv1);
860 tlist2.AddToList(&fill3);
861 tlist2.AddToList(&fill4);
862 tlist2.AddToList(&fill2a);
863 tlist2.AddToList(&fill3a);
864 tlist2.AddToList(&fill4a);
865 tlist2.AddToList(&fill5a);
866 tlist2.AddToList(&fill6a);
867 tlist2.AddToList(&fill7a);
868 tlist2.AddToList(&fill8a);
869 //tlist2.AddToList(&fill9a);
870
871 MEvtLoop loop2(fName);
872 loop2.SetParList(&plist);
873 loop2.SetDisplay(fDisplay);
874 loop2.SetLogStream(fLog);
875
876 if (!SetupEnv(loop2))
877 return kFALSE;
878
879 if (!loop2.Eventloop(fMaxEvents))
880 {
881 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
882 return kFALSE;
883 }
884
885 if (!loop2.GetDisplay())
886 {
887 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
888 return kFALSE;
889 }
890
891 // -------------------------- Spectrum ----------------------------
892
893 DisplaySpectrum(area, excess, hest, ontime);
894 DisplaySize(plist);
895
896 if (!fPathOut.IsNull())
897 fDisplay->SaveAsRoot(fPathOut);
898
899 return kTRUE;
900}
Note: See TracBrowser for help on using the repository browser.