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

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