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

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