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

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