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

Last change on this file since 9337 was 9335, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 57.5 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-2008
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 <TLine.h>
39#include <TFile.h>
40#include <TGraph.h>
41#include <TChain.h>
42#include <TLatex.h>
43#include <TCanvas.h>
44#include <TObjArray.h>
45
46// Environment
47#include "MLog.h"
48#include "MLogManip.h"
49
50#include "MMath.h"
51#include "MString.h"
52#include "MDirIter.h"
53
54#include "MStatusArray.h"
55#include "MStatusDisplay.h"
56
57// Container
58#include "MH3.h"
59#include "MHn.h"
60#include "MBinning.h"
61#include "MParameters.h"
62#include "MDataSet.h"
63#include "MMcCorsikaRunHeader.h"
64
65// Spectrum
66#include "../mhflux/MAlphaFitter.h"
67#include "../mhflux/MHAlpha.h"
68#include "../mhflux/MHCollectionArea.h"
69#include "../mhflux/MHEnergyEst.h"
70#include "../mhflux/MMcSpectrumWeight.h"
71
72// Eventloop
73#include "MEvtLoop.h"
74#include "MTaskList.h"
75#include "MParList.h"
76
77// Tasks/Filter
78#include "MReadMarsFile.h"
79#include "MEnergyEstimate.h"
80#include "MTaskEnv.h"
81#include "MFillH.h"
82#include "MFDataPhrase.h"
83#include "MFDataMember.h"
84#include "MContinue.h"
85#include "MWriteRootFile.h"
86
87ClassImp(MJSpectrum);
88
89using namespace std;
90
91MJSpectrum::MJSpectrum(const char *name, const char *title)
92 : fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0),
93 fEstimateEnergy(0), fCalcHadronness(0), fCalcDisp(0), fForceTheta(kFALSE)
94{
95 fName = name ? name : "MJSpectrum";
96 fTitle = title ? title : "Standard program to calculate spectrum";
97
98 // Make sure that fDisplay is maintained properly
99 // (i.e. removed via RecursiveRemove if closed)
100 gROOT->GetListOfCleanups()->Add(this);
101 SetBit(kMustCleanup);
102}
103
104// --------------------------------------------------------------------------
105//
106// Delete all the fCut* data members and fCalcHadronness/fEstimateEnergy
107//
108MJSpectrum::~MJSpectrum()
109{
110 if (fCutQ)
111 delete fCutQ;
112 if (fCut0)
113 delete fCut0;
114 if (fCut1)
115 delete fCut1;
116 if (fCut2)
117 delete fCut2;
118 if (fCut3)
119 delete fCut3;
120 if (fCutS)
121 delete fCutS;
122 if (fEstimateEnergy)
123 delete fEstimateEnergy;
124 if (fCalcHadronness)
125 delete fCalcHadronness;
126 if (fCalcDisp)
127 delete fCalcDisp;
128}
129
130// --------------------------------------------------------------------------
131//
132// Setup a task estimating the energy. The given task is cloned.
133//
134void MJSpectrum::SetEnergyEstimator(const MTask *task)
135{
136 if (fEstimateEnergy)
137 delete fEstimateEnergy;
138 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
139}
140
141// --------------------------------------------------------------------------
142//
143// Read a MTask named name into task from the open file. If this task is
144// required mustexist can be set. Depending on success kTRUE or kFALSE is
145// returned. If the task is a MContinue SetAllowEmpty is called.
146// The name of the task is set to name.
147//
148Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name, Bool_t mustexist) const
149{
150 // If a task is set delete task
151 if (task)
152 {
153 delete task;
154 task = 0;
155 }
156
157 // Try to get task from file
158 task = (MTask*)gFile->Get(name);
159 if (!task)
160 {
161 if (!mustexist)
162 return kTRUE;
163 *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
164 return kFALSE;
165 }
166
167 // Check if task inherits from MTask
168 if (!task->InheritsFrom(MTask::Class()))
169 {
170 *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
171 delete task;
172 return kFALSE;
173 }
174
175 // Set name of task
176 task->SetName(name);
177
178 // SetAllowEmpty for MContinue tasks
179 if (dynamic_cast<MContinue*>(task))
180 dynamic_cast<MContinue*>(task)->SetAllowEmpty();
181
182 return kTRUE;
183}
184
185// --------------------------------------------------------------------------
186//
187// Print setup used for the MC processing, namely MAlphaFitter ans all fCut*.
188//
189void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
190{
191 fLog->Separator("Alpha Fitter");
192 *fLog << all;
193 fit.Print("result");
194
195 fLog->Separator("Used Cuts");
196 fCutQ->Print();
197 fCut0->Print();
198 fCut1->Print();
199 fCutS->Print();
200 fCut2->Print();
201 fCut3->Print();
202}
203
204// --------------------------------------------------------------------------
205//
206// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
207// the dataset.
208// The simulated energy range and spectral slope is initialized from
209// there.
210// In the following eventloops the forced check in MMcSpectrumWeight
211// ensures, that the spectral slope and energy range doesn't change.
212//
213Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
214{
215 fLog->Separator("Initialize energy weighting");
216
217 // Read Resources
218 if (!CheckEnv(w))
219 {
220 *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
221 return kFALSE;
222 }
223
224 w.Print("new");
225
226 return kTRUE;
227}
228
229Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
230{
231 *fLog << inf << "Reading from file: " << fPathIn << endl;
232
233 TFile file(fPathIn, "READ");
234 if (!file.IsOpen())
235 {
236 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
237 return -1;
238 }
239
240 MStatusArray arr;
241 if (arr.Read()<=0)
242 {
243 *fLog << err << dbginf << "ERROR - MStatusDisplay not found in file... abort." << endl;
244 return -1;
245 }
246
247 //arr.EnableTH1Workaround();
248
249 TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
250 TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
251 if (!vstime)
252 {
253 *fLog << err << dbginf << "ERROR - Theta [TH1D] not found in OnTime-tab... abort." << endl;
254 return -1;
255 }
256 if (!size)
257 {
258 *fLog << err << dbginf << "ERROR - Excess [TH1D] not found in Hist-tab... abort." << endl;
259 return -1;
260 }
261
262 vstime->Copy(h1);
263 size->Copy(h2);
264 h1.SetDirectory(0);
265 h2.SetDirectory(0);
266
267 if (fDisplay)
268 arr.DisplayIn(*fDisplay, "Hist");
269
270 if (!ReadTask(fCutQ, "CutQ"))
271 return -1;
272 if (!ReadTask(fCut0, "Cut0"))
273 return -1;
274 if (!ReadTask(fCut1, "Cut1"))
275 return -1;
276 if (!ReadTask(fCut2, "Cut2"))
277 return -1;
278 if (!ReadTask(fCut3, "Cut3"))
279 return -1;
280 if (!ReadTask(fCalcHadronness, "CalcHadronness", kFALSE))
281 return -1;
282 if (!ReadTask(fCalcDisp, "CalcDisp", kFALSE))
283 return -1;
284
285 TObjArray arrread;
286
287 TIter Next(plist);
288 TObject *o=0;
289 while ((o=Next()))
290 if (o->InheritsFrom(MBinning::Class()))
291 arrread.Add(o);
292
293 arrread.Add(plist.FindObject("MAlphaFitter"));
294
295 if (!ReadContainer(arrread))
296 return -1;
297
298 const Double_t ufl = vstime->GetBinContent(0);
299 const Double_t ofl = vstime->GetBinContent(vstime->GetNbinsX()+1);
300 const Double_t eff = vstime->Integral()+ufl+ofl;
301 if (ufl>0)
302 {
303 if (vstime->GetBinLowEdge(1)<=0)
304 {
305 *fLog << err << "ERROR - Undeflow bin of OnTime histogram not empty as it ought to be." << endl;
306 return -1;
307 }
308 *fLog << warn << "WARNING - " << Form("%.1f%% (%.0fs)", 100*ufl/eff, ufl) << " of the eff. observation time is in the underflow bin." << endl;
309 }
310 if (ofl>0)
311 *fLog << warn << "WARNING - " << Form("%.1f%% (%.0fs)", 100*ofl/eff, ofl) << " of the eff. observation time is in the overflow bin." << endl;
312
313 if (!fForceRunTime)
314 return eff;
315
316 TH1D *time = (TH1D*)arr.FindObjectInCanvas("ExcessTime", "TH1D", "Hist");
317 if (!time)
318 {
319 *fLog << err;
320 *fLog << "ERROR - ExcessTime [TH1D] not found in Hist-tab... abort." << endl;
321 *fLog << " Did you try to process Monte Carlos with --force-runtime?" <<endl;
322 return -1;
323 }
324
325 const Double_t obs = time->GetXaxis()->GetXmax()-time->GetXaxis()->GetXmin();
326
327 *fLog << inf;
328 *fLog << "Total run time: " << obs/60 << "min" << endl;
329 *fLog << "Eff. obs. time: " << eff/60 << "min (" << Form("%.1f%%", 100*eff/obs) << ")" << endl;
330
331 return obs;
332}
333
334// --------------------------------------------------------------------------
335//
336// return maximum value of MMcRunHeader.fImpactMax stored in the RunHeaders
337// of the files from the MC dataset
338//
339Bool_t MJSpectrum::AnalyzeMC(const MDataSet &set, Float_t &impactmax, Float_t &emin/*, Float_t emax*/) const
340{
341 if (fDisplay)
342 fDisplay->SetStatusLine1("Analyzing Monte Carlo headers...");
343
344 // ----- Create chain -----
345 *fLog << inf << "Getting files... " << flush;
346 TChain chain("RunHeaders");
347 if (!set.AddFilesOn(chain))
348 return kFALSE;
349 *fLog << "done. " << endl;
350
351 *fLog << all;
352 *fLog << "Searching for maximum impact... " << flush;
353 impactmax = chain.GetMaximum("MMcRunHeader.fImpactMax");
354 *fLog << "found " << impactmax/100 << "m" << endl;
355
356 *fLog << "Searching for minimum lower energy limit... " << flush;
357 emin = chain.GetMinimum("MMcCorsikaRunHeader.fELowLim");
358 *fLog << "found " << emin << "GeV" << endl;
359
360 *fLog << "Searching for maximum lower energy limit... " << flush;
361 const Float_t emin2 = chain.GetMaximum("MMcCorsikaRunHeader.fELowLim");
362 *fLog << "found " << emin2 << "GeV" << endl;
363
364 // Need a check for the upper energy LIMIT?!?
365
366 if (emin2>emin)
367 {
368 *fLog << warn;
369 *fLog << "WARNING - You are using files with different lower limits for the simulated" << endl;
370 *fLog << " energy, i.e. that the collection area and your correction" << endl;
371 *fLog << " coefficients between " << emin << "GeV and ";
372 *fLog << emin2 << "GeV might be wrong." << endl;
373 }
374
375/*
376 *fLog << "Getting maximum energy... " << flush;
377 emax = chain.GetMaximum("MMcCorsikaRunHeader.fEUppLim");
378 *fLog << "found " << emax << "GeV" << endl;
379 */
380 return kTRUE;
381}
382
383Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
384{
385 // Some debug output
386 *fLog << all << endl;
387 fLog->Separator("Compiling original MC distribution");
388
389 Float_t impactmax=0, Emin=0;
390 if (!AnalyzeMC(set, impactmax, Emin))
391 return kFALSE;
392
393 *fLog << inf << "Getting files... " << flush;
394 MDirIter iter;
395 if (!set.AddFilesOn(iter))
396 return kFALSE;
397 *fLog << "done. " << endl;
398
399 const Int_t tot = iter.GetNumEntries();
400
401 // Prepare histogram
402 h.Reset();
403 h.Sumw2();
404 if (h.InheritsFrom(TH2::Class()))
405 {
406 h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
407 h.SetXTitle("\\Theta [\\circ]");
408 h.SetYTitle("E [GeV]");
409 h.SetZTitle(MString::Format("Counts normalized to r=%.0fm", impactmax/100));
410 }
411 else
412 {
413 h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
414 h.SetXTitle("\\Theta [\\circ]");
415 h.SetYTitle(MString::Format("Counts normalized to r=%.0fm", impactmax/100));
416 }
417 h.SetDirectory(0);
418
419 const TString cont = h.InheritsFrom(TH2::Class()) ?
420 "MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC" :
421 "MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC";
422
423 if (fDisplay)
424 fDisplay->SetStatusLine1("Reading OriginalMC distribution...");
425
426 TH1 *hfill = (TH1*)h.Clone(h.InheritsFrom(TH2::Class())?"ThetaEMC":"ThetaMC");
427 hfill->SetDirectory(0);
428
429 *fLog << all << "Compiling simulated source spectrum... please stand by, this may take a while." << endl;
430
431 Double_t evts = 0;
432 Int_t num = 0;
433
434 // Reading this with a eventloop is five times slower :(
435 TString fname;
436 while (fDisplay)
437 {
438 if (fDisplay)
439 fDisplay->SetProgressBarPosition(Float_t(num++)/tot);
440
441 // Get next filename
442 fname = iter.Next();
443 if (fname.IsNull())
444 break;
445
446 if (fDisplay)
447 fDisplay->SetStatusLine2(fname);
448
449 // open file
450 TFile file(fname);
451 if (file.IsZombie())
452 {
453 *fLog << err << "ERROR - Couldn't open file " << fname << endl;
454 return kFALSE;
455 }
456
457 // Get tree OriginalMC
458 TTree *t = dynamic_cast<TTree*>(file.Get("OriginalMC"));
459 if (!t)
460 {
461 *fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl;
462 return kFALSE;
463 }
464 if (t->GetEntries()==0)
465 {
466 *fLog << warn << "WARNING - OriginalMC of " << fname << " empty." << endl;
467 continue;
468 }
469
470 // Get tree RunHeaders
471 TTree *rh = dynamic_cast<TTree*>(file.Get("RunHeaders"));
472 if (!rh)
473 {
474 *fLog << err << "ERROR - File " << fname << " doesn't contain tree RunHeaders." << endl;
475 return kFALSE;
476 }
477 if (rh->GetEntries()!=1)
478 {
479 *fLog << err << "ERROR - RunHeaders of " << fname << " doesn't contain exactly one entry." << endl;
480 return kFALSE;
481 }
482
483 // Get corsika run header
484 MMcCorsikaRunHeader *head=0;
485 rh->SetBranchAddress("MMcCorsikaRunHeader.", &head);
486 rh->GetEntry(0);
487 if (!head)
488 {
489 *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from " << fname << "." << endl;
490 return kFALSE;
491 }
492
493 // Get the maximum impact parameter of this file. Due to different
494 // production areas an additional scale-factor is applied.
495 //
496 // Because it is assumed that the efficiency outside the production
497 // area is nearly zero no additional weight has to be applied to the
498 // events after cuts. For the events before cuts it is fair to use
499 // weights... maybe filling the residual impact with unweighted
500 // events would be better?!? (Not that the weighting might be
501 // less correct with low statistics, because it could pronounce
502 // a fluctuation)
503 const Double_t impact = rh->GetMaximum("MMcRunHeader.fImpactMax");
504 const Double_t scale = impactmax/impact;
505
506 // Propagate the run header to MMcSpectrumWeight
507 if (!weight.Set(*head))
508 {
509 *fLog << err << "ERROR - Initializing MMcSpectrumWeight from " << fname << " failed." << endl;
510 return kFALSE;
511 }
512
513 // Get the corresponding rule for the weights
514 const TString weights(weight.GetFormulaWeights("MMcEvtBasic"));
515
516 // No we found everything... go on reading contents
517 *fLog << inf2 << "Reading OriginalMC of " << fname << endl;
518
519 // Fill histogram from tree
520 hfill->SetDirectory(&file);
521 if (t->Draw(cont, weights, "goff")<0)
522 {
523 *fLog << err << "ERROR - Reading OriginalMC failed." << endl;
524 return kFALSE;
525 }
526 hfill->SetDirectory(0);
527
528 evts += hfill->GetEntries();
529
530 // ----- Complete incomplete energy ranges -----
531 // ----- and apply production area weights -----
532 weight.CompleteEnergySpectrum(*hfill, Emin, scale*scale);
533
534 // Add weighted data from file to final histogram
535 h.Add(hfill);
536 }
537 delete hfill;
538
539 *fLog << all << "Total number of events from files in Monte Carlo dataset: " << evts << endl;
540 if (fDisplay)
541 fDisplay->SetStatusLine2("done.");
542
543 if (!fDisplay || h.GetEntries()>0)
544 return kTRUE;
545
546 *fLog << err << "ERROR - Histogram with distribution from OriginalMC empty..." << endl;
547 return kFALSE;
548}
549
550void MJSpectrum::GetThetaDistribution(TH1D &temp1, TH2D &hist2) const
551{
552 TH1D &temp2 = *hist2.ProjectionX();
553
554 // Display some stuff
555 if (fDisplay)
556 {
557 TCanvas &c = fDisplay->AddTab("ZdDist");
558 c.Divide(2,2);
559
560 // On-Time vs. Theta
561 c.cd(1);
562 gPad->SetBorderMode(0);
563 gPad->SetGridx();
564 gPad->SetGridy();
565 temp1.DrawCopy();
566
567 // Number of MC events (produced) vs Theta
568 c.cd(2);
569 gPad->SetBorderMode(0);
570 gPad->SetGridx();
571 gPad->SetGridy();
572 temp2.SetName("NVsTheta");
573 temp2.DrawCopy("hist");
574
575 c.cd(4);
576 gPad->SetBorderMode(0);
577 gPad->SetGridx();
578 gPad->SetGridy();
579
580 c.cd(3);
581 gPad->SetBorderMode(0);
582 gPad->SetGridx();
583 gPad->SetGridy();
584 }
585
586 // Calculate the Probability
587 temp1.Divide(&temp2);
588 temp1.Scale(1./temp1.Integral(1, temp1.GetNbinsX()+1));
589
590 // Some cosmetics: Name, Axis, etc.
591 temp1.SetName("ProbVsTheta");
592 temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
593 temp1.SetYTitle("Probability");
594 if (fDisplay)
595 temp1.DrawCopy("hist");
596
597 delete &temp2;
598}
599
600// --------------------------------------------------------------------------
601//
602// Display the final theta distribution.
603//
604Bool_t MJSpectrum::DisplayResult(const TH2D &h2) const
605{
606 if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
607 {
608 *fLog << err << "ERROR - Display or tab ZdDist vanished... abort." << endl;
609 return kFALSE;
610 }
611
612 TH1D &proj = *h2.ProjectionX();
613 proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
614 proj.SetXTitle("\\Theta [\\circ]");
615 proj.SetYTitle("Counts");
616 proj.SetLineColor(kBlue);
617 proj.SetDirectory(0);
618 proj.SetBit(kCanDelete);
619
620 TVirtualPad *pad = gPad;
621
622 pad->cd(4);
623 proj.DrawCopy();
624
625 pad->cd(1);
626 TH1D *theta = (TH1D*)gPad->FindObject("Theta");
627 if (!theta)
628 {
629 *fLog << err << "ERROR - Theta-Histogram vanished... cannot proceed." << endl;
630 return kFALSE;
631 }
632
633 // Check whether histogram is empty
634 if (proj.GetMaximum()==0)
635 {
636 *fLog << err;
637 *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap" << endl;
638 *fLog << " with the Zenith Angle distribution of your observation." << endl;
639 *fLog << " Maybe the energy binning is undefined or wrong (from ";
640 *fLog << h2.GetYaxis()->GetXmin() << "GeV to " << h2.GetYaxis()->GetXmax() << "GeV)" << endl;
641 theta->SetLineColor(kRed);
642 return kFALSE;;
643 }
644
645 // scale histogram and set new maximum for display
646 proj.Scale(theta->GetMaximum()/proj.GetMaximum());
647 theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
648
649 // draw project
650 proj.Draw("same");
651
652 // Compare both histograms
653 *fLog << inf << "Comparing theta distributions for data and MCs." << endl;
654
655 const Double_t prob = proj.Chi2Test(theta, "P");
656 if (prob==1)
657 return kTRUE;
658
659 if (prob>0.99)
660 {
661 *fLog << inf;
662 *fLog << "The Zenith Angle distribution of your Monte Carlos fits well" << endl;
663 *fLog << "with the Zenith Angle distribution of your observation." << endl;
664 *fLog << "The probability for identical Theta distributions is " << prob << endl;
665 return kTRUE;
666 }
667
668 if (prob<0.01)
669 {
670 *fLog << err;
671 *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit" << endl;
672 *fLog << " with the Zenith Angle distribution of your observation." << endl;
673 *fLog << " The probability for identical Theta distributions is " << prob << endl;
674 if (!fForceTheta)
675 *fLog << " To force processing use --force-theta (with care!)" << endl;
676 theta->SetLineColor(kRed);
677 return fForceTheta;
678 }
679
680 *fLog << warn;
681 *fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well" << endl;
682 *fLog << " with the Zenith Angle distribution of your observation." << endl;
683 *fLog << " The probability for identical Theta distributions is " << prob << endl;
684 return kTRUE;
685}
686
687TString MJSpectrum::GetHAlpha() const
688{
689 TString cls("MHAlpha");
690 if (!fDisplay)
691 return cls;
692
693 TCanvas *c = fDisplay->GetCanvas("Hist");
694 if (!c)
695 return cls;
696
697 TIter Next(c->GetListOfPrimitives());
698 TObject *obj=0;
699 while ((obj=Next()))
700 if (obj->InheritsFrom(MHAlpha::Class()))
701 break;
702
703 return obj ? TString(obj->ClassName()) : cls;
704}
705
706// --------------------------------------------------------------------------
707//
708// Fills the excess histogram (vs E-est) from the events stored in the
709// ganymed result file and therefor estimates the energy.
710//
711// The resulting histogram excess-vs-energy ist copied into h2.
712//
713Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2)/*const*/
714{
715 // Now fill the histogram
716 *fLog << endl;
717 fLog->Separator("Refill Excess");
718 *fLog << endl;
719
720 MTaskList tlist;
721 plist.AddToList(&tlist);
722
723 MReadTree read("Events");
724 read.DisableAutoScheme();
725 read.AddFile(fPathIn);
726/*
727 MTaskEnv taskenv0("CalcDisp");
728 taskenv0.SetDefault(fCalcDisp);
729
730 MTaskEnv taskenv1("CalcHadronness");
731 taskenv1.SetDefault(fCalcHadronness);
732 */
733 MEnergyEstimate est;
734 MTaskEnv taskenv2("EstimateEnergy");
735 taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
736
737 MContinue *cont = new MContinue("", "CutS");
738 cont->SetAllowEmpty();
739
740 if (fCutS)
741 delete fCutS;
742 fCutS = cont;
743
744 // Try to find the class used to determine the signal!
745 const TString cls = GetHAlpha();
746
747 // FIXME: Create HistE and HistEOff to be able to modify it from
748 // the resource file.
749 MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
750 MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
751
752 MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
753 MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
754
755 fill1.SetFilter(&f0);
756 fill2.SetFilter(&f1);
757
758 tlist.AddToList(&read);
759 //tlist.AddToList(&taskenv0); // not necessary, stored in file!
760 //tlist.AddToList(&taskenv1); // not necessary, stored in file!
761 tlist.AddToList(fCutS);
762 tlist.AddToList(&taskenv2);
763 tlist.AddToList(&f0);
764 tlist.AddToList(&f1);
765 tlist.AddToList(&fill1);
766 tlist.AddToList(&fill2);
767
768 // by setting it here it is distributed to all consecutive tasks
769 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
770
771 MEvtLoop loop("RefillExcess"); // ***** fName *****
772 loop.SetParList(&plist);
773 loop.SetDisplay(fDisplay);
774 loop.SetLogStream(fLog);
775
776 if (!SetupEnv(loop))
777 return kFALSE;
778
779 if (!loop.Eventloop())
780 {
781 *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
782 return kFALSE;
783 }
784
785 if (!loop.GetDisplay())
786 {
787 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
788 return kFALSE;
789 }
790
791 const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE");
792 if (!halpha)
793 {
794 *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
795 return kFALSE;
796 }
797
798 halpha->GetHEnergy().Copy(h2);
799 h2.SetDirectory(0);
800
801 return kTRUE;
802}
803
804/*
805Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
806{
807 MTaskList tlist1;
808 plist.Replace(&tlist1);
809
810 MReadMarsFile readmc("OriginalMC");
811 //readmc.DisableAutoScheme();
812 if (!set.AddFilesOn(readmc))
813 return kFALSE;
814
815 readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
816 readmc.EnableBranch("MMcEvtBasic.fEnergy");
817
818 mh1.SetLogy();
819 mh1.SetLogz();
820 mh1.SetName("ThetaE");
821
822 MFillH fill0(&mh1);
823 //fill0.SetDrawOption("projx only");
824
825 MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
826 MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
827 if (bins2 && bins3)
828 {
829 bins2->SetName("BinningThetaEY");
830 bins3->SetName("BinningThetaEX");
831 }
832 tlist1.AddToList(&readmc);
833 tlist1.AddToList(&weight);
834
835 temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
836 MH3 mh3mc(temp1);
837
838 MFEventSelector2 sel1(mh3mc);
839 sel1.SetHistIsProbability();
840
841 fill0.SetFilter(&sel1);
842
843 //if (!fRawMc)
844 tlist1.AddToList(&sel1);
845 tlist1.AddToList(&fill0);
846
847 // by setting it here it is distributed to all consecutive tasks
848 tlist1.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
849
850 MEvtLoop loop1("IntermediateLoop"); // ***** fName *****
851 loop1.SetParList(&plist);
852 loop1.SetLogStream(fLog);
853 loop1.SetDisplay(fDisplay);
854
855 if (!SetupEnv(loop1))
856 return kFALSE;
857
858 if (!loop1.Eventloop(fMaxEvents))
859 {
860 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
861 return kFALSE;
862 }
863
864 if (!loop1.GetDisplay())
865 {
866 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
867 return kFALSE;
868 }
869
870 if (bins2 && bins3)
871 {
872 bins2->SetName("BinningEnergyEst");
873 bins3->SetName("BinningTheta");
874 }
875
876 return kTRUE;
877}
878*/
879
880TString MJSpectrum::FormFloat(Double_t d)
881{
882 TString s;
883 s += d;
884 return s.Strip(TString::kLeading);
885}
886
887TString MJSpectrum::FormFlux(const TF1 &f, const char *unit)
888{
889 Double_t p0 = -f.GetParameter(0);
890 Double_t p1 = f.GetParameter(1);
891
892 Double_t e0 = f.GetParError(0);
893 Double_t e1 = f.GetParError(1);
894
895 MMath::Format(p0, e0);
896 MMath::Format(p1, e1);
897
898 const Int_t i = TMath::FloorNint(TMath::Log10(p1));
899 const Double_t exp = TMath::Power(10, i);
900
901 TString str = MString::Format("(%s #pm %s)·10^{%d} ",
902 FormFloat(p1/exp).Data(), FormFloat(e1/exp).Data(), i);
903
904 str += MString::Format("#left(#frac{E}{%s}#right)^{-%s #pm %s}", unit,
905 FormFloat(p0).Data(), FormFloat(e0).Data());
906
907 str += " TeV^{-1} m^{-2} s^{-1}";
908
909 return str;
910}
911
912TString MJSpectrum::FormString(const TF1 &f, Byte_t type)
913{
914 switch (type)
915 {
916 case 0:
917 return FormFlux(f, "500GeV");
918 case 1:
919 return MString::Format("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF());
920 case 2:
921 return MString::Format("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(), f.GetNDF()));
922 }
923 return "";
924}
925
926TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
927{
928 Axis_t lo, hi;
929 MH::GetRangeUser(spectrum, lo, hi);
930
931 TF1 f("f", "[1]*(x/500)^[0]", lo, hi);
932 f.SetParameter(0, -3.0);
933 f.SetParameter(1, spectrum.GetMaximum());
934 f.SetLineColor(kBlue);
935 f.SetLineWidth(2);
936 spectrum.Fit(&f, "NIR"); // M skipped
937 f.DrawCopy("same");
938
939 TString str = FormString(f);
940
941 TLatex tex;
942 tex.SetTextSize(0.045);
943 tex.SetBit(TLatex::kTextNDC);
944 tex.SetTextAlign(31);
945 tex.DrawLatex(0.89, 0.935, str);
946
947 str = FormString(f, 1);
948 tex.DrawLatex(0.89, 0.83, str);
949
950 str = FormString(f, 2);
951 tex.DrawLatex(0.89, 0.735, str);
952
953 TArrayD res(2);
954 res[0] = f.GetParameter(0);
955 res[1] = f.GetParameter(1);
956 return res;
957}
958
959// --------------------------------------------------------------------------
960//
961// Calculate the final spectrum from:
962// - collection area
963// - excess
964// - correction coefficients
965// - ontime
966// and display it
967//
968TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
969{
970 // Create copies of the histograms
971 TH1D collarea(area.GetHEnergy());
972 TH1D spectrum(excess);
973 TH1D weights;
974
975 // Get spill-over corrections from energy estimation
976 hest.GetWeights(weights); // E_mc/E_est
977
978 // Print effective on-time
979 cout << "Effective On time: " << ontime << "s" << endl;
980
981 // Setup spectrum plot
982 spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
983 spectrum.SetYTitle("N/sm^{2}");
984 spectrum.SetDirectory(NULL);
985 spectrum.SetBit(kCanDelete);
986
987 // Divide by collection are and on-time
988 spectrum.Scale(1./ontime);
989 spectrum.Divide(&collarea);
990
991 // Draw spectrum before applying spill-over corrections
992 TCanvas &c1 = fDisplay->AddTab("Spectrum");
993 c1.Divide(2,2);
994 c1.cd(1);
995 gPad->SetBorderMode(0);
996 gPad->SetLogx();
997 gPad->SetLogy();
998 gPad->SetGridx();
999 gPad->SetGridy();
1000 collarea.DrawCopy();
1001
1002 c1.cd(2);
1003 gPad->SetBorderMode(0);
1004 gPad->SetLogx();
1005 gPad->SetLogy();
1006 gPad->SetGridx();
1007 gPad->SetGridy();
1008 TH1D *spec=(TH1D*)spectrum.DrawCopy();
1009 //FitSpectrum(*spec);
1010
1011 c1.cd(3);
1012 gPad->SetBorderMode(0);
1013 gPad->SetLogx();
1014 gPad->SetLogy();
1015 gPad->SetGridx();
1016 gPad->SetGridy();
1017 weights.DrawCopy();
1018
1019 // Apply spill-over correction (done't take the errors into account)
1020 // They are supposed to be identical with the errors of the
1021 // collection area and cancel out.
1022 //spectrum.Multiply(&weights);
1023 spectrum.SetNameTitle("Flux", "Spectrum");
1024 spectrum.SetBit(TH1::kNoStats);
1025
1026 // Minimum number of excessevents to get 3sigma in 1h
1027 TF1 sensl("SensLZA", "85*(x/200)^(-0.55)", 100, 6000);
1028 TF1 sensh("SensHZA", "35*(x/200)^(-0.70)", 100, 1000);
1029 //sens.SetLineColor(kBlue);
1030 //sens.DrawClone("Csame");
1031
1032 c1.cd(4);
1033 gPad->SetBorderMode(0);
1034 gPad->SetLogx();
1035 gPad->SetLogy();
1036 gPad->SetGridx();
1037 gPad->SetGridy();
1038
1039 TGraph gsensl;//("Sensitivity LZA", "");
1040 TGraph gsensh;//("Sensitivity HZA", "");
1041
1042 const Double_t sqrton = TMath::Sqrt(ontime/3600.);
1043
1044 for (int i=0; i<excess.GetNbinsX(); i++)
1045 {
1046 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
1047 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1));
1048
1049 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
1050 spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
1051
1052 if (collarea.GetBinContent(i+1)<=0)
1053 continue;
1054
1055 const Double_t cen = spectrum.GetBinCenter(i+1);
1056 gsensl.SetPoint(gsensl.GetN(), cen, sensl.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime);
1057 gsensh.SetPoint(gsensh.GetN(), cen, sensh.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime);
1058
1059 cout << cen << " " << sensl.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime;
1060 cout << " " << sensh.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime;
1061 cout << endl;
1062 }
1063
1064 spectrum.SetMinimum(1e-12);
1065 spectrum.SetXTitle("E [GeV]");
1066 spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
1067 spec = (TH1D*)spectrum.DrawCopy();
1068
1069 TLatex tex;
1070 tex.SetTextColor(13);
1071
1072 TF1 fc("Crab", "7.0e-6*(x/300)^(-2.31-0.26*log10(x/300))", 100, 6000);
1073 fc.SetLineStyle(9);
1074 fc.SetLineWidth(2);
1075 fc.SetLineColor(14);
1076 fc.DrawCopy("same");
1077
1078 tex.DrawLatex(1250, fc.Eval(1250), "Crab/\\Gamma=-2.3");
1079
1080 TF1 fa("PG1553", "1.8e-6*(x/200)^-4.21", 90, 600);
1081 static_cast<const TAttLine&>(fc).Copy(fa);
1082 fa.DrawCopy("same");
1083
1084 tex.SetTextAlign(23);
1085 tex.DrawLatex(600, fa.Eval(600), "PG1553/\\Gamma=-4.2");
1086
1087 gROOT->SetSelectedPad(0);
1088
1089 gsensl.SetLineStyle(5);
1090 gsensl.SetLineColor(14);
1091 gsensl.SetLineWidth(2);
1092 gsensl.DrawClone("C")->SetBit(kCanDelete);
1093
1094 gsensh.SetLineStyle(5);
1095 gsensh.SetLineColor(14);
1096 gsensh.SetLineWidth(2);
1097 gsensh.DrawClone("C")->SetBit(kCanDelete);
1098
1099 // Display dN/dE*E^2 for conveinience
1100 fDisplay->AddTab("Check");
1101 gPad->SetLogx();
1102 gPad->SetLogy();
1103 gPad->SetGrid();
1104
1105 // Calculate Spectrum * E^2
1106 for (int i=0; i<spectrum.GetNbinsX(); i++)
1107 {
1108 const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)*spectrum.GetBinLowEdge(i+2))*1e-3;
1109
1110 spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e);
1111 spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) *e*e);
1112 }
1113
1114 for (int i=0; i<gsensl.GetN(); i++)
1115 {
1116 const Double_t e = gsensl.GetX()[i]*1e-3;
1117
1118 gsensl.GetY()[i] *= e*e;
1119 gsensh.GetY()[i] *= e*e;
1120 }
1121
1122 spectrum.SetName("FluxStd");
1123 spectrum.SetMarkerStyle(kFullDotMedium);
1124 spectrum.SetTitle("Differential flux times E^{2}");
1125 spectrum.SetYTitle("E^{2}·dN/dE [N·TeV/sm^{2}]");
1126 spectrum.SetDirectory(0);
1127 spectrum.DrawCopy();
1128
1129 TF1 fc2("Crab*E^2", "x^2*Crab*1e-6", 100, 6000);
1130 static_cast<const TAttLine&>(fc).Copy(fc2);
1131 fc2.DrawCopy("same");
1132
1133 TF1 fa2("PG1553*E^2", "x^2*PG1553*1e-6", 100, 6000);
1134 static_cast<const TAttLine&>(fc).Copy(fa2);
1135 fa2.DrawCopy("same");
1136
1137 gsensl.DrawClone("C")->SetBit(kCanDelete);
1138 gsensh.DrawClone("C")->SetBit(kCanDelete);
1139
1140 // Fit Spectrum
1141 c1.cd(4);
1142 return FitSpectrum(*spec);
1143
1144/*
1145 TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
1146 f.SetParameter(0, -2.87);
1147 f.SetParameter(1, 1.9e-6);
1148 f.SetLineColor(kGreen);
1149 spectrum.Fit(&f, "NIM", "", 100, 5000);
1150 f.DrawCopy("same");
1151
1152 const Double_t p0 = f.GetParameter(0);
1153 const Double_t p1 = f.GetParameter(1);
1154
1155 const Double_t e0 = f.GetParError(0);
1156 const Double_t e1 = f.GetParError(1);
1157
1158 const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
1159 const Double_t exp = TMath::Power(10, np);
1160
1161 TString str;
1162 str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
1163 str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
1164 str += "\\frac{ph}{TeVm^{2}s}";
1165
1166 TLatex tex;
1167 tex.SetTextSize(0.045);
1168 tex.SetBit(TLatex::kTextNDC);
1169 tex.DrawLatex(0.45, 0.935, str);
1170
1171 str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
1172 tex.DrawLatex(0.70, 0.83, str);
1173
1174 TArrayD res(2);
1175 res[0] = f.GetParameter(0);
1176 res[1] = f.GetParameter(1);
1177
1178 return res;
1179 */
1180}
1181
1182// --------------------------------------------------------------------------
1183//
1184// Scale some image parameter plots using the scale factor and plot them
1185// together with the corresponding MC histograms.
1186// Called from DisplaySize
1187//
1188Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
1189{
1190 TString same(name);
1191 same += "Same";
1192
1193 TH1 *h1 = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
1194 TH1 *h2 = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
1195 if (!h1 || !h2)
1196 return kFALSE;
1197
1198 TObject *obj = plist.FindObject(plot);
1199 if (!obj)
1200 {
1201 *fLog << warn << plot << " not in parameter list... skipping." << endl;
1202 return kFALSE;
1203 }
1204
1205 TH1 *h3 = (TH1*)obj->FindObject(name);
1206 if (!h3)
1207 {
1208 *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
1209 return kFALSE;
1210 }
1211
1212
1213 const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
1214 const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
1215
1216 gPad->SetBorderMode(0);
1217 h2->SetLineColor(kBlack);
1218 h3->SetLineColor(kBlue);
1219 h2->Add(h1, -ascale);
1220
1221 //h2->Scale(1./ontime); //h2->Integral());
1222 h3->Scale(scale); //h3->Integral());
1223
1224 h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
1225
1226 h2 = h2->DrawCopy();
1227 h3 = h3->DrawCopy("same");
1228
1229 // Don't do this on the original object!
1230 h2->SetStats(kFALSE);
1231 h3->SetStats(kFALSE);
1232
1233 return kTRUE;
1234}
1235
1236// --------------------------------------------------------------------------
1237//
1238// Take a lot of histograms and plot them together in one plot.
1239// Calls PlotSame
1240//
1241Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
1242{
1243 *fLog << inf << "Reading from file: " << fPathIn << endl;
1244
1245 TFile file(fPathIn, "READ");
1246 if (!file.IsOpen())
1247 {
1248 *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
1249 return kFALSE;
1250 }
1251
1252 file.cd();
1253 MStatusArray arr;
1254 if (arr.Read()<=0)
1255 {
1256 *fLog << "MStatusDisplay not found in file... abort." << endl;
1257 return kFALSE;
1258 }
1259
1260 TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
1261 if (!excess)
1262 return kFALSE;
1263
1264 // ------------------- Plot excess versus size -------------------
1265
1266 TCanvas &c = fDisplay->AddTab("Excess");
1267 c.Divide(3,2);
1268 c.cd(1);
1269 gPad->SetBorderMode(0);
1270 gPad->SetLogx();
1271 gPad->SetLogy();
1272 gPad->SetGridx();
1273 gPad->SetGridy();
1274
1275 excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
1276 excess = excess->DrawCopy("E2");
1277 // Don't do this on the original object!
1278 excess->SetStats(kFALSE);
1279 excess->SetMarkerStyle(kFullDotMedium);
1280 excess->SetFillColor(kBlack);
1281 excess->SetFillStyle(0);
1282 excess->SetName("Excess ");
1283 excess->SetDirectory(0);
1284
1285 TObject *o=0;
1286 if ((o=plist.FindObject("ExcessMC")))
1287 {
1288 TH1 *histsel = (TH1F*)o->FindObject("");
1289 if (histsel)
1290 {
1291 if (scale<0)
1292 scale = excess->Integral()/histsel->Integral();
1293
1294 histsel->Scale(scale);
1295 histsel->SetLineColor(kBlue);
1296 histsel->SetBit(kCanDelete);
1297 histsel = histsel->DrawCopy("E1 same");
1298 // Don't do this on the original object!
1299 histsel->SetStats(kFALSE);
1300 /*
1301 fLog->Separator("Kolmogorov Test");
1302 histsel->KolmogorovTest(excess, "DX");
1303 fLog->Separator("Chi^2 Test");
1304 const Double_t p = histsel->Chi2Test(excess, "P");
1305
1306 TLatex tex;
1307 tex.SetBit(TLatex::kTextNDC);
1308 tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
1309 */
1310 }
1311 }
1312
1313 // -------------- Comparison of Image Parameters --------------
1314 c.cd(2);
1315 PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale);
1316
1317 c.cd(3);
1318 PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
1319
1320 c.cd(4);
1321 PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale);
1322
1323 c.cd(5);
1324 PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale);
1325
1326 c.cd(6);
1327 PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale);
1328
1329 return kTRUE;
1330}
1331
1332// --------------------------------------------------------------------------
1333//
1334void MJSpectrum::DisplayCutEfficiency(const MHCollectionArea &area0, const MHCollectionArea &area1) const
1335{
1336 if (!fDisplay)
1337 return;
1338
1339 const TH1D &trig = area0.GetHEnergy();
1340 TH1D &cut = (TH1D&)*area1.GetHEnergy().Clone();
1341
1342 fDisplay->AddTab("CutEff");
1343
1344 gPad->SetBorderMode(0);
1345 gPad->SetFrameBorderMode(0);
1346 gPad->SetLogx();
1347 gPad->SetGridx();
1348 gPad->SetGridy();
1349
1350 cut.Divide(&trig);
1351 cut.Scale(100);
1352 cut.SetNameTitle("CutEff", "Background Supression: Cut efficiency (after star)");
1353 cut.SetYTitle("\\eta [%]");
1354 cut.SetDirectory(0);
1355 cut.SetMinimum(0);
1356 cut.SetMaximum(100);
1357 cut.SetBit(kCanDelete);
1358 cut.Draw();
1359
1360 TLine line;
1361 line.SetLineColor(kBlue);
1362 line.SetLineWidth(2);
1363 line.SetLineStyle(kDashed);
1364 line.DrawLine(cut.GetBinLowEdge(1), 50, cut.GetBinLowEdge(cut.GetNbinsX()+1), 50);
1365}
1366
1367void MJSpectrum::SetupHistEvtDist(MHn &hist) const
1368{
1369 hist.AddHist("MMcEvt.fEnergy");
1370 hist.InitName("EnergyDist;EnergyEst");
1371 hist.InitTitle("Unweighted event distribution (Real statistics);E [GeV];Counts;");
1372
1373 hist.AddHist("MPointingPos.fZd");
1374 hist.InitName("ZdDist;Theta");
1375 hist.InitTitle("Unweighted event distribution (Real statistics);Zd [\\circ];Counts;");
1376}
1377
1378void MJSpectrum::SetupHistEnergyEst(MHn &hist) const
1379{
1380 const char *res = "log10(MEnergyEst.fVal)-log10(MMcEvt.fEnergy)";
1381
1382 hist.AddHist("MHillas.fSize", res);
1383 hist.InitName("ResSize;Size;EnergyResidual");
1384 hist.InitTitle(";S [phe];\\Delta lg E;");
1385 hist.SetDrawOption("colz profx");
1386
1387 hist.AddHist("MPointingPos.fZd", res);
1388 hist.InitName("ResTheta;Theta;EnergyResidual");
1389 hist.InitTitle(";Zd [\\circ];\\Delta lg E;");
1390 hist.SetDrawOption("colz profx");
1391
1392 hist.AddHist("MNewImagePar.fLeakage1", res);
1393 hist.InitName("ResLeak;Leakage;EnergyResidual");
1394 hist.InitTitle(";Leak;\\Delta lg E;");
1395 hist.SetDrawOption("colz profx");
1396
1397 hist.AddHist("MHillasSrc.fDist*3.37e-3", res);
1398 hist.InitName("ResDist;Dist;EnergyResidual");
1399 hist.InitTitle(";D [\\circ];\\Delta lg E;");
1400 hist.SetDrawOption("colz profx");
1401}
1402
1403void MJSpectrum::SetupHistDisp(MHn &hist) const
1404{
1405 const char *res = "-Disp.fVal*sign(MHillasSrc.fCosDeltaAlpha)-MHillasSrc.fDist*3.37e-3";
1406
1407 hist.AddHist("MHillas.fSize", res);
1408 hist.InitName("ResSize;Size;ResidualDist");
1409 hist.InitTitle(";S [phe];Disp-Dist [\\circ];");
1410 hist.SetDrawOption("colz profx");
1411
1412 hist.AddHist("MPointingPos.fZd", res);
1413 hist.InitName("ResTheta;Theta;ResidualDist");
1414 hist.InitTitle(";Zd [\\circ];Disp-Dist [\\circ];");
1415 hist.SetDrawOption("colz profx");
1416
1417 hist.AddHist("MNewImagePar.fLeakage1", res);
1418 hist.InitName("ResLeak;Leakage;ResidualDist");
1419 hist.InitTitle(";Leak;Disp-Dist [\\circ];");
1420 hist.SetDrawOption("colz profx");
1421
1422 hist.AddHist("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/3.37e-3", res);
1423 hist.InitName("ResSlope;Slope;ResidualDist");
1424 hist.InitTitle(";Slope;Disp-Dist [\\circ];");
1425 hist.SetDrawOption("colz profx");
1426}
1427
1428void MJSpectrum::SetupHistEnergyRes(MHn &hist) const
1429{
1430 hist.AddHist("MEnergyEst.fVal", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
1431 hist.InitName("ResEest;EnergyEst;");
1432 hist.InitTitle(";E_{est} [GeV];Resolution (E_{mc}/E_{est}-1)^{2};");
1433
1434 //hist.AddHist("MMcEvt.fEnergy", "(MEnergyEst.fVal/MMcEvt.fEnergy-1)^2", MH3::kProfile);
1435 //hist.InitName("ResEmc;EnergyEst;");
1436 //hist.InitTitle(";E_{mc} [GeV];Resolution (E_{est}/E_{mc}-1)^{2};");
1437 hist.AddHist("MHillas.fSize", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
1438 hist.InitName("ResSize;Size;");
1439 hist.InitTitle(";S [phe];Resolution (E_{mc}/E_{est}-1)^{2};");
1440
1441 hist.AddHist("MPointingPos.fZd", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
1442 hist.InitName("ResTheta;Theta;");
1443 hist.InitTitle(";\\Theta [\\circ];Resolution (E_{mc}/E_{est}-1)^{2};");
1444
1445 hist.AddHist("MMcEvt.fImpact/100", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
1446 hist.InitName("ResImpact;Impact;");
1447 hist.InitTitle(";I [m];Resolution (E_{mc}/E_{est}-1)^{2};");
1448}
1449
1450// --------------------------------------------------------------------------
1451//
1452// Setup write to write:
1453// container tree optional?
1454// -------------- ---------- -----------
1455// "MHillas" to "Events"
1456// "MHillasSrc" to "Events"
1457// "Hadronness" to "Events" yes
1458// "MEnergyEst" to "Events" yes
1459// "DataType" to "Events"
1460//
1461void MJSpectrum::SetupWriter(MWriteRootFile *write/*, const char *name*/) const
1462{
1463 if (!write)
1464 return;
1465
1466 //write->SetName(name);
1467 write->AddContainer("MHillas", "Events");
1468 write->AddContainer("MHillasSrc", "Events");
1469 write->AddContainer("MHillasExt", "Events");
1470 //write->AddContainer("MPointingPos", "Events");
1471 write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
1472 write->AddContainer("MImagePar", "Events", kFALSE);
1473 write->AddContainer("MNewImagePar", "Events", kFALSE);
1474 write->AddContainer("MNewImagePar2", "Events", kFALSE);
1475 write->AddContainer("Hadronness", "Events", kFALSE);
1476 write->AddContainer("MSrcPosCam", "Events", kFALSE);
1477 write->AddContainer("MSrcPosAnti", "Events", kFALSE);
1478 write->AddContainer("ThetaSquared", "Events", kFALSE);
1479 write->AddContainer("OpticalAxis", "Events", kFALSE);
1480 write->AddContainer("Disp", "Events", kFALSE);
1481 write->AddContainer("Ghostbuster", "Events", kFALSE);
1482 write->AddContainer("MEnergyEst", "Events", kFALSE);
1483 write->AddContainer("MTime", "Events", kFALSE);
1484 write->AddContainer("MMcEvt", "Events", kFALSE);
1485 write->AddContainer("MWeight", "Events");
1486 write->AddContainer("DataType", "Events");
1487 write->AddContainer("FileId", "Events");
1488 write->AddContainer("EvtNumber", "Events");
1489}
1490
1491Bool_t MJSpectrum::Process(const MDataSet &set)
1492{
1493 if (!set.IsValid())
1494 {
1495 *fLog << err << "ERROR - DataSet invalid!" << endl;
1496 return kFALSE;
1497 }
1498
1499 if (!HasWritePermission(GetPathOut()))
1500 return kFALSE;
1501
1502 if (!CheckEnv())
1503 return kFALSE;
1504
1505 // --------------------------------------------------------------------------------
1506
1507 *fLog << inf;
1508 fLog->Separator(GetDescriptor());
1509 *fLog << "Compile Monte Carlo sample (dataset " << set.GetBaseName() << ")" << endl;
1510 *fLog << endl;
1511
1512 if (fDisplay)
1513 fDisplay->SetWindowName(fName);
1514
1515 // Setup everything which is read from the ganymed file
1516 MBinning bins1("BinningAlpha");
1517 MBinning bins2("BinningEnergyEst");
1518 MBinning bins3("BinningTheta");
1519 MBinning bins4("BinningFalseSource");
1520 MBinning bins5("BinningWidth");
1521 MBinning bins6("BinningLength");
1522 MBinning bins7("BinningDist");
1523 MBinning bins8("BinningM3Long");
1524 MBinning bins9("BinningM3Trans");
1525 MBinning bins0("BinningSlope");
1526 MBinning binsa("BinningAsym");
1527 MBinning binsb("BinningConc1");
1528
1529 MAlphaFitter fit;
1530
1531 MParList plist;
1532 plist.AddToList(&bins0);
1533 plist.AddToList(&bins1);
1534 plist.AddToList(&bins3);
1535 plist.AddToList(&bins4);
1536 plist.AddToList(&bins5);
1537 plist.AddToList(&bins6);
1538 plist.AddToList(&bins7);
1539 plist.AddToList(&bins8);
1540 plist.AddToList(&bins9);
1541 plist.AddToList(&binsa);
1542 plist.AddToList(&binsb);
1543 plist.AddToList(&fit);
1544
1545 // Read from the ganymed file
1546 TH1D htheta, size;
1547 const Float_t ontime = ReadInput(plist, htheta, size);
1548 if (ontime<0)
1549 {
1550 *fLog << err << GetDescriptor() << ": Could not determine effective on time..." << endl;
1551 return kFALSE;
1552 }
1553
1554 // Set Zenith angle binning to binning from the ganymed-file
1555 bins3.SetEdges(htheta, 'x');
1556
1557 // Read energy binning from resource file
1558 if (!CheckEnv(bins2))
1559 {
1560 *fLog << err << "ERROR - Reading energy binning from resources failed." << endl;
1561 return kFALSE;
1562 }
1563 plist.AddToList(&bins2); // For later use in MC processing
1564
1565 // Initialize weighting to a new spectrum as defined in the resource file
1566 MMcSpectrumWeight weight;
1567 if (!InitWeighting(set, weight))
1568 return kFALSE;
1569
1570 // Print Theta and energy binning for cross-checks
1571 *fLog << all << endl;
1572 bins2.Print();
1573 bins3.Print();
1574
1575 // Now we read the MC distribution as produced by corsika
1576 // vs zenith angle and energy.
1577 // Weight for the new energy spectrum defined in MMcSpectumWeight
1578 // are applied.
1579 // Also correction for different lower energy bounds and
1580 // different production areas (impact parameters) are applied.
1581 TH2D hist;
1582 hist.UseCurrentStyle();
1583 MH::SetBinning(&hist, &bins3, &bins2);
1584 if (!ReadOrigMCDistribution(set, hist, weight))
1585 return kFALSE;
1586
1587 // Check if user has closed the display
1588 if (!fDisplay)
1589 return kTRUE;
1590
1591 // Display histograms and calculate za-weights into htheta
1592 GetThetaDistribution(htheta, hist);
1593
1594 // Give the zenith angle weights to the weighting task
1595 weight.SetWeightsZd(&htheta);
1596
1597 // No we apply the the zenith-angle-weights to the corsika produced
1598 // MC distribution. Unfortunately this must be done manually
1599 // because we are multiplying column by column
1600 for (int y=0; y<=hist.GetNbinsY()+1; y++)
1601 for (int x=0; x<=hist.GetNbinsX()+1; x++)
1602 {
1603 hist.SetBinContent(x, y, hist.GetBinContent(x, y)*htheta.GetBinContent(x));
1604 hist.SetBinError(x, y, hist.GetBinError(x, y) *htheta.GetBinContent(x));
1605 }
1606
1607 // Display the resulting distribution and check it matches
1608 // the observation time distribution (this could be false
1609 // for example if you miss MCs of some zenith angles, which you have
1610 // data for)
1611 if (!DisplayResult(hist))
1612 return kFALSE;
1613
1614 // -------------- Fill excess events versus energy ---------------
1615
1616 // Refill excess histogram to determine the excess events
1617 TH1D excess;
1618 if (!Refill(plist, excess))
1619 return kFALSE;
1620
1621 // Print the setup and result of the MAlphaFitter, print used cuts
1622 PrintSetup(fit);
1623
1624
1625 // ------------------------- Final loop --------------------------
1626
1627 *fLog << endl;
1628 fLog->Separator("Calculate efficiencies");
1629 *fLog << endl;
1630
1631 MTaskList tlist2;
1632 plist.AddToList(&tlist2);
1633
1634 MReadMarsFile read("Events");
1635 read.DisableAutoScheme();
1636 if (!set.AddFilesOn(read))
1637 return kFALSE;
1638
1639 // Selector to get correct (final) theta-distribution
1640 //temp1.SetXTitle("MPointingPos.fZd");
1641 //
1642 //MH3 mh3(temp1);
1643 //
1644 //MFEventSelector2 sel2(mh3);
1645 //sel2.SetHistIsProbability();
1646 //
1647 //MContinue contsel(&sel2);
1648 //contsel.SetInverted();
1649
1650 // Get correct source position
1651 //MSrcPosCalc calc;
1652
1653 // Calculate corresponding Hillas parameters
1654 //MHillasCalc hcalc1;
1655 //MHillasCalc hcalc2("MHillasCalcAnti");
1656 //hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
1657 //hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
1658 //hcalc2.SetNameHillasSrc("MHillasSrcAnti");
1659 //hcalc2.SetNameSrcPosCam("MSrcPosAnti");
1660
1661 // Fill collection area and energy estimator (unfolding)
1662 // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
1663 MHCollectionArea area0("TriggerArea");
1664 MHCollectionArea area1;
1665 area0.SetHistAll(hist);
1666 area1.SetHistAll(hist);
1667
1668 MHEnergyEst hest;
1669
1670 MFillH fill30(&area0, "", "FillTriggerArea");
1671 MFillH fill31(&area1, "", "FillCollectionArea");
1672 MFillH fill4(&hest, "", "FillEnergyEst");
1673 MFillH fill5("MHThreshold", "", "FillThreshold");
1674 fill30.SetWeight();
1675 fill31.SetWeight();
1676 fill4.SetWeight();
1677 fill5.SetWeight();
1678 fill30.SetNameTab("TrigArea");
1679 fill31.SetNameTab("ColArea");
1680 fill4.SetNameTab("E-Est");
1681 fill5.SetNameTab("Threshold");
1682
1683/*
1684 MH3 henergy("MMcEvt.fEnergy");
1685 henergy.SetName("EventDist;EnergyEst");
1686 henergy.SetTitle("Unweighted event distribution (Real statistics);E [GeV];Counts;");
1687 henergy.Sumw2();
1688 */
1689
1690 // ---------------------------------------------------------
1691
1692 MBinning binsA(50, 10, 100000, "BinningSize", "log");
1693 MBinning binsC(50, 0, 0.3, "BinningLeakage", "lin");
1694 MBinning binsB(51, -1, 1, "BinningEnergyResidual", "lin");
1695 MBinning binsD(51, -1, 1, "BinningResidualDist", "lin");
1696 MBinning binsI(16, 0, 800, "BinningImpact", "lin");
1697
1698 plist.AddToList(&binsA);
1699 plist.AddToList(&binsB);
1700 plist.AddToList(&binsC);
1701 plist.AddToList(&binsD);
1702 plist.AddToList(&binsI);
1703
1704 MHn heest("Energy", "Energy Residual (lg E_{est} - lg E_{mc})");
1705 SetupHistEnergyEst(heest);
1706
1707 MHn hdisp("Disp", "Dist residual (Disp-Dist)");
1708 SetupHistDisp(hdisp);
1709
1710 MHn henergy("EvtDist");
1711 SetupHistEvtDist(henergy);
1712
1713 MHn heres("EnergyRes");
1714 SetupHistEnergyRes(heres);
1715
1716 MFillH fill4b(&heest, "", "FillEnergyResidual");
1717 fill4b.SetWeight();
1718
1719 MFillH fill4c(&hdisp, "", "FillDispResidual");
1720 fill4c.SetWeight();
1721
1722 MFillH fill4d(&heres, "", "FillEnergyResolution");
1723 fill4d.SetWeight();
1724
1725 MFDataPhrase fdisp("Disp.fVal*sign(MHillasSrc.fCosDeltaAlpha)<0", "FilterDisp");
1726 fill4c.SetFilter(&fdisp);
1727
1728 // ---------------------------------------------------------
1729
1730 MH3 hsize("MHillas.fSize");
1731 hsize.SetName("ExcessMC");
1732 hsize.Sumw2();
1733
1734 MBinning bins(size, "BinningExcessMC");
1735 plist.AddToList(&hsize);
1736 plist.AddToList(&bins);
1737
1738 // ---------------------------------------------------------
1739
1740 MFillH fillsp("MHSrcPosCam", "MSrcPosCam", "FillSrcPosCam");
1741 MFillH fill0a(&henergy, "", "FillEventDist");
1742 MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
1743 MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
1744 MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
1745 MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
1746 MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
1747 MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
1748 MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
1749 MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC");
1750 fillsp.SetNameTab("SrcPos");
1751 fill0a.SetNameTab("EvtDist");
1752 fill1a.SetNameTab("PreCut");
1753 fill2a.SetNameTab("PostCut");
1754 fill3a.SetNameTab("VsSize");
1755 fill4a.SetNameTab("HilExt");
1756 fill5a.SetNameTab("HilSrc");
1757 fill6a.SetNameTab("ImgPar");
1758 fill7a.SetNameTab("NewPar");
1759 fill8a.SetBit(MFillH::kDoNotDisplay);
1760 fill1a.SetWeight();
1761 fill2a.SetWeight();
1762 fill3a.SetWeight();
1763 fill4a.SetWeight();
1764 fill5a.SetWeight();
1765 fill6a.SetWeight();
1766 fill7a.SetWeight();
1767 fill8a.SetWeight();
1768 fillsp.SetWeight();
1769
1770 // Try to find the class used to determine the signal!
1771 // FIXME: To be done: A task checking the lower 1% after the lower
1772 // energy limit!
1773
1774 MTaskEnv taskenv0("CalcDisp");
1775 taskenv0.SetDefault(fCalcDisp);
1776
1777 MTaskEnv taskenv1("CalcHadronness");
1778 taskenv1.SetDefault(fCalcHadronness);
1779
1780 MEnergyEstimate est;
1781 MTaskEnv taskenv2("EstimateEnergy");
1782 taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
1783
1784 MWriteRootFile write(GetPathOut());
1785 SetupWriter(&write);
1786
1787 MParameterI *par = (MParameterI*)plist.FindCreateObj("MParameterI", "DataType");
1788 if (!par)
1789 return kFALSE;
1790 par->SetVal(2);
1791
1792 tlist2.AddToList(&read);
1793 // If no weighting should be applied but the events should
1794 // be thrown away according to the theta distribution
1795 // it is enabled here
1796 //if (!fRawMc && fNoThetaWeights)
1797 // tlist2.AddToList(&contsel);
1798 //tlist2.AddToList(&calc);
1799 //tlist2.AddToList(&hcalc1);
1800 //tlist2.AddToList(&hcalc2);
1801 tlist2.AddToList(&weight);
1802 tlist2.AddToList(&fillsp);
1803 tlist2.AddToList(&fill1a);
1804 tlist2.AddToList(&fill30);
1805 tlist2.AddToList(fCutQ);
1806 tlist2.AddToList(fCut0);
1807 tlist2.AddToList(&taskenv0);
1808 tlist2.AddToList(&taskenv1);
1809 tlist2.AddToList(&fdisp);
1810 tlist2.AddToList(&fill4c);
1811 tlist2.AddToList(fCut1);
1812 tlist2.AddToList(fCutS);
1813 tlist2.AddToList(fCut2);
1814 tlist2.AddToList(fCut3);
1815 tlist2.AddToList(&taskenv2);
1816 if (!GetPathOut().IsNull())
1817 tlist2.AddToList(&write);
1818 tlist2.AddToList(&fill31);
1819 tlist2.AddToList(&fill4);
1820 tlist2.AddToList(&fill4b);
1821 tlist2.AddToList(&fill4d);
1822 tlist2.AddToList(&fill5);
1823 tlist2.AddToList(&fill0a);
1824 tlist2.AddToList(&fill2a);
1825 tlist2.AddToList(&fill3a);
1826 tlist2.AddToList(&fill4a);
1827 tlist2.AddToList(&fill5a);
1828 tlist2.AddToList(&fill6a);
1829 tlist2.AddToList(&fill7a);
1830 tlist2.AddToList(&fill8a);
1831 //tlist2.AddToList(&fill9a);
1832
1833 // by setting it here it is distributed to all consecutive tasks
1834 tlist2.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
1835
1836 MEvtLoop loop2(fName); // ***** fName *****
1837 loop2.SetParList(&plist);
1838 loop2.SetDisplay(fDisplay);
1839 loop2.SetLogStream(fLog);
1840
1841 if (!SetupEnv(loop2))
1842 return kFALSE;
1843
1844 if (!loop2.Eventloop(fMaxEvents))
1845 {
1846 *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
1847 return kFALSE;
1848 }
1849
1850 if (!loop2.GetDisplay())
1851 {
1852 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
1853 return kFALSE;
1854 }
1855
1856 fLog->Separator("Energy Estimator");
1857 if (plist.FindObject("EstimateEnergy"))
1858 plist.FindObject("EstimateEnergy")->Print();
1859
1860 fLog->Separator("Spectrum");
1861
1862 // -------------------------- Spectrum ----------------------------
1863
1864 // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
1865 TArrayD res(DisplaySpectrum(area1, excess, hest, ontime));
1866
1867 // Spectrum fitted (convert res[1] from TeV to GeV)
1868 TF1 flx("flx", MString::Format("%e*pow(x/500, %f)", res[1]/500, res[0]));
1869
1870 // Number of events this spectrum would produce per s and m^2
1871 Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
1872
1873 // scale with effective collection area to get the event rate (N/s)
1874 // scale with the effective observation time to absolute observed events
1875 n *= area1.GetCollectionAreaAbs()*ontime; // N
1876
1877 // Now calculate the scale factor from the number of events
1878 // produced and the number of events which should have been
1879 // observed with our telescope in the time ontime
1880 const Double_t scale = n/area1.GetEntries();
1881
1882 // Print normalization constant
1883 cout << "MC normalization factor: " << scale << endl;
1884
1885 // Display cut efficiency
1886 DisplayCutEfficiency(area0, area1);
1887
1888 // Overlay normalized plots
1889 DisplaySize(plist, scale);
1890
1891 // check if output should be written
1892 if (!fPathOut.IsNull())
1893 {
1894 TNamed ganame("ganymed.root", fPathIn.Data());
1895
1896 // Write the output
1897 TObjArray cont;
1898 cont.Add((TObject*)GetEnv()); // const_cast
1899 cont.Add((TObject*)&set); // const_cast
1900 cont.Add(plist.FindObject("MAlphaFitter"));
1901 cont.Add(&area0);
1902 cont.Add(&area1);
1903 cont.Add(&hest);
1904 cont.Add(&ganame);
1905
1906 if (fDisplay)
1907 cont.Add(fDisplay);
1908
1909 if (!WriteContainer(cont, "", GetPathOut().IsNull()?"RECREATE":"UPDATE"))
1910 {
1911 *fLog << err << GetDescriptor() << ": Writing result failed." << endl;
1912 return kFALSE;
1913 }
1914 }
1915
1916 *fLog << all << GetDescriptor() << ": Done." << endl;
1917 *fLog << endl << endl;
1918
1919 return kTRUE;
1920}
Note: See TracBrowser for help on using the repository browser.