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

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