Context Navigation

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

Visit:

Last change on this file since 8646 was 8587, checked in by tbretz, 17 years ago
* empty log message *
File size: 37.2 KB

Line
1	/* ======================================================================== *\
2	!
3	! *
4	! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5	! * Software. It is distributed to you in the hope that it can be a useful
6	! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7	! * It is distributed WITHOUT ANY WARRANTY.
8	! *
9	! * Permission to use, copy, modify and distribute this software and its
10	! * documentation for any purpose is hereby granted without fee,
11	! * provided that the above copyright notice appear in all copies and
12	! * that both that copyright notice and this permission notice appear
13	! * in supporting documentation. It is provided "as is" without express
14	! * or implied warranty.
15	! *
16	!
17	!
18	! Author(s): Thomas Bretz, 4/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19	!
20	! Copyright: MAGIC Software Development, 2000-2005
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 "MStatusArray.h"
50	#include "MStatusDisplay.h"
51
52	// Container
53	#include "MH3.h"
54	#include "MBinning.h"
55	#include "MDataSet.h"
56	#include "MMcCorsikaRunHeader.h"
57
58	// Spectrum
59	#include "../mhflux/MAlphaFitter.h"
60	#include "../mhflux/MHAlpha.h"
61	#include "../mhflux/MHCollectionArea.h"
62	//#include "../mhflux/MHThreshold.h"
63	#include "../mhflux/MHEnergyEst.h"
64	#include "../mhflux/MMcSpectrumWeight.h"
65
66	// Eventloop
67	#include "MEvtLoop.h"
68	#include "MTaskList.h"
69	#include "MParList.h"
70
71	// Tasks/Filter
72	#include "MReadMarsFile.h"
73	#include "MReadMarsFile.h"
74	#include "MFEventSelector2.h"
75	#include "MFDataMember.h"
76	#include "MEnergyEstimate.h"
77	#include "MTaskEnv.h"
78	#include "MFillH.h"
79	#include "MHillasCalc.h"
80	//#include "MSrcPosCalc.h"
81	#include "MContinue.h"
82
83	ClassImp(MJSpectrum);
84
85	using namespace std;
86
87	MJSpectrum::MJSpectrum(const char name, const char title)
88	: fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
89	fCalcHadronness(0), fRefill(kFALSE), fSimpleMode(kTRUE),
90	fRawMc(kFALSE), fNoThetaWeights(kFALSE)
91	{
92	fName = name ? name : "MJSpectrum";
93	fTitle = title ? title : "Standard program to calculate spectrum";
94	}
95
96	MJSpectrum::~MJSpectrum()
97	{
98	if (fCut0)
99	delete fCut0;
100	if (fCut1)
101	delete fCut1;
102	if (fCut2)
103	delete fCut2;
104	if (fCut3)
105	delete fCut3;
106	if (fEstimateEnergy)
107	delete fEstimateEnergy;
108	if (fCalcHadronness)
109	delete fCalcHadronness;
110	}
111
112	// --------------------------------------------------------------------------
113	//
114	// Setup a task estimating the energy. The given task is cloned.
115	//
116	void MJSpectrum::SetEnergyEstimator(const MTask *task)
117	{
118	if (fEstimateEnergy)
119	delete fEstimateEnergy;
120	fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
121	}
122
123	Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name, Bool_t mustexist) const
124	{
125	if (task)
126	{
127	delete task;
128	task = 0;
129	}
130
131	task = (MTask*)gFile->Get(name);
132	if (!task)
133	{
134	if (!mustexist)
135	return kTRUE;
136	*fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
137	return kFALSE;
138	}
139	if (!task->InheritsFrom(MTask::Class()))
140	{
141	*fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
142	delete task;
143	return kFALSE;
144	}
145
146	task->SetName(name);
147
148	if (dynamic_cast<MContinue*>(task))
149	dynamic_cast<MContinue*>(task)->SetAllowEmpty();
150
151	return kTRUE;
152	}
153
154	void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
155	{
156	fLog->Separator("Alpha Fitter");
157	*fLog << all;
158	fit.Print();
159
160	fLog->Separator("Used Cuts");
161	fCut0->Print();
162	fCut1->Print();
163	fCut2->Print();
164	fCut3->Print();
165
166	//gLog.Separator("Energy Estimator");
167	//fEstimateEnergy->Print();
168	}
169
170	// --------------------------------------------------------------------------
171	//
172	// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
173	// the dataset.
174	// The simulated energy range and spectral slope is initialized from
175	// there.
176	// In the following eventloops the forced check in MMcSpectrumWeight
177	// ensures, that the spectral slope and energy range doesn't change.
178	//
179	Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
180	{
181	fLog->Separator("Initialize energy weighting");
182
183	// Read Resources
184	if (!CheckEnv(w))
185	{
186	*fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
187	return kFALSE;
188	}
189
190	// Read the first corsika RunHeader from the MC files
191	TChain chain("RunHeaders");
192	if (!set.AddFilesOn(chain))
193	return kFALSE;
194
195	MMcCorsikaRunHeader *h=0;
196	chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
197	chain.GetEntry(1);
198
199	if (!h)
200	{
201	*fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
202	return kFALSE;
203	}
204
205	// Propagate the run header to MMcSpectrumWeight
206	if (!w.Set(*h))
207	{
208	*fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
209	return kFALSE;
210	}
211
212	// Print new setup of MMcSpectrumWeight
213	w.Print();
214	return kTRUE;
215	}
216
217	Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
218	{
219	*fLog << inf << "Reading from file: " << fPathIn << endl;
220
221	TFile file(fPathIn, "READ");
222	if (!file.IsOpen())
223	{
224	*fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
225	return -1;
226	}
227
228	MStatusArray arr;
229	if (arr.Read()<=0)
230	{
231	*fLog << "MStatusDisplay not found in file... abort." << endl;
232	return -1;
233	}
234
235	TH1D vstime = (TH1D)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
236	TH1D size = (TH1D)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
237	if (!vstime \|\| !size)
238	return -1;
239
240	vstime->Copy(h1);
241	size->Copy(h2);
242	h1.SetDirectory(0);
243	h2.SetDirectory(0);
244
245	if (fDisplay)
246	arr.DisplayIn(*fDisplay, "Hist");
247
248	if (!ReadTask(fCut0, "Cut0"))
249	return -1;
250	if (!ReadTask(fCut1, "Cut1"))
251	return -1;
252	if (!ReadTask(fCut2, "Cut2"))
253	return -1;
254	if (!ReadTask(fCut3, "Cut3"))
255	return -1;
256	if (!ReadTask(fCalcHadronness, "CalcHadronness", kFALSE))
257	return -1;
258
259	TObjArray arrread;
260
261	TIter Next(plist);
262	TObject *o=0;
263	while ((o=Next()))
264	if (o->InheritsFrom(MBinning::Class()))
265	arrread.Add(o);
266
267	arrread.Add(plist.FindObject("MAlphaFitter"));
268
269	if (!ReadContainer(arrread))
270	return -1;
271
272	return vstime->Integral();
273	}
274
275	Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
276	{
277	// Some debug output
278	fLog->Separator("Compiling original MC distribution");
279
280	// The name of the input container is MMcEvtBasic
281	weight.SetNameMcEvt("MMcEvtBasic");
282	// Get the corresponding rule for the weights
283	const TString w(weight.GetFormulaWeights());
284	// Set back to MMcEvt
285	weight.SetNameMcEvt();
286
287	*fLog << inf << "Using weights: " << w << endl;
288	*fLog << "Please stand by, this may take a while..." << flush;
289
290	if (fDisplay)
291	fDisplay->SetStatusLine1("Compiling MC distribution...");
292
293	// Create chain
294	TChain chain("OriginalMC");
295	if (!set.AddFilesOn(chain))
296	return kFALSE;
297
298	// Prepare histogram
299	h.Reset();
300	h.Sumw2();
301
302	// Fill histogram from chain
303	h.SetDirectory(gROOT);
304	if (h.InheritsFrom(TH2::Class()))
305	{
306	h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
307	h.SetXTitle("\\Theta [\\circ]");
308	h.SetYTitle("E [GeV]");
309	h.SetZTitle("Counts");
310	chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
311	}
312	else
313	{
314	h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
315	h.SetXTitle("\\Theta [\\circ]");
316	h.SetYTitle("Counts");
317	chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
318	}
319	h.SetDirectory(0);
320
321	*fLog << "done." << endl;
322	if (fDisplay)
323	fDisplay->SetStatusLine2("done.");
324
325	if (h.GetEntries()>0)
326	return kTRUE;
327
328	*fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
329
330	return h.GetEntries()>0;
331	}
332
333	Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
334	{
335	// Display some stuff
336	if (fDisplay)
337	{
338	TCanvas &c = fDisplay->AddTab("ZdDist");
339	c.Divide(2,2);
340
341	// On-Time vs. Theta
342	c.cd(1);
343	gPad->SetBorderMode(0);
344	temp1.DrawCopy();
345
346	// Number of MC events (produced) vs Theta
347	c.cd(2);
348	gPad->SetBorderMode(0);
349	temp2.SetName("NVsTheta");
350	temp2.DrawCopy("hist");
351
352	c.cd(4);
353	gPad->SetBorderMode(0);
354
355	c.cd(3);
356	gPad->SetBorderMode(0);
357	}
358
359	// Calculate the Probability
360	temp1.Divide(&temp2);
361	temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral());
362
363	// Some cosmetics: Name, Axis, etc.
364	temp1.SetName("ProbVsTheta");
365	temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
366	temp1.SetYTitle("Probability");
367	if (fDisplay)
368	temp1.DrawCopy("hist");
369
370	return kTRUE;
371	}
372
373	// --------------------------------------------------------------------------
374	//
375	// Display the final theta distribution.
376	//
377	void MJSpectrum::DisplayResult(const TH2D &h2) const
378	{
379	if (!fDisplay \|\| !fDisplay->CdCanvas("ZdDist"))
380	return;
381
382	TH1D &proj = *h2.ProjectionX();
383	proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
384	proj.SetXTitle("\\Theta [\\circ]");
385	proj.SetYTitle("Counts");
386	proj.SetLineColor(kBlue);
387	proj.SetDirectory(0);
388	proj.SetBit(kCanDelete);
389
390	TVirtualPad *pad = gPad;
391
392	pad->cd(4);
393	proj.DrawCopy();
394
395	pad->cd(1);
396	TH1D theta = (TH1D)gPad->FindObject("Theta");
397	if (theta)
398	{
399	proj.Scale(theta->GetMaximum()/proj.GetMaximum());
400	theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
401	}
402	proj.Draw("same");
403	}
404
405	// --------------------------------------------------------------------------
406	//
407	// Fills the excess histogram (vs E-est) from the events stored in the
408	// ganymed result file and therefor estimates the energy.
409	//
410	// The resulting histogram excess-vs-energy ist copied into h2.
411	//
412	Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
413	{
414	// Try to find the class used to determin the signal!
415	TString cls("MHAlpha");
416	if (fDisplay)
417	{
418	TCanvas *c = fDisplay->GetCanvas("Hist");
419	if (c)
420	{
421	TIter Next(c->GetListOfPrimitives());
422	TObject *obj=0;
423	while ((obj=Next()))
424	if (obj->InheritsFrom(MHAlpha::Class()))
425	break;
426	if (obj)
427	cls = obj->ClassName();
428	}
429	}
430
431	cout << "FOUND: "<< cls << endl;
432
433	// Now fill the histogram
434	*fLog << endl;
435	fLog->Separator("Refill Excess");
436	*fLog << endl;
437
438	MTaskList tlist;
439	plist.AddToList(&tlist);
440
441	MReadTree read("Events");
442	read.DisableAutoScheme();
443	read.AddFile(fPathIn);
444
445	MEnergyEstimate est;
446	MTaskEnv taskenv1("EstimateEnergy");
447	taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
448
449	// FIXME: Create HistE and HistEOff to be able to modify it from
450	// the resource file.
451
452	MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
453	MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
454
455	MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
456	MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
457
458	fill1.SetFilter(&f0);
459	fill2.SetFilter(&f1);
460
461	tlist.AddToList(&read);
462	tlist.AddToList(&taskenv1);
463	tlist.AddToList(&f0);
464	tlist.AddToList(&f1);
465	tlist.AddToList(&fill1);
466	tlist.AddToList(&fill2);
467
468	// by setting it here it is distributed to all consecutive tasks
469	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
470
471	MEvtLoop loop("RefillExcess"); // *** fName ***
472	loop.SetParList(&plist);
473	loop.SetDisplay(fDisplay);
474	loop.SetLogStream(fLog);
475
476	if (!SetupEnv(loop))
477	return kFALSE;
478
479	if (!loop.Eventloop())
480	{
481	*fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
482	return kFALSE;
483	}
484
485	if (!loop.GetDisplay())
486	{
487	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
488	return kFALSE;
489	}
490
491	const MHAlpha halpha = (MHAlpha )plist.FindObject("HistE");
492	if (!halpha)
493	{
494	*fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
495	return kFALSE;
496	}
497
498	halpha->GetHEnergy().Copy(h2);
499	h2.SetDirectory(0);
500
501	return kTRUE;
502	}
503
504	Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
505	{
506	MTaskList tlist1;
507	plist.Replace(&tlist1);
508
509	MReadMarsFile readmc("OriginalMC");
510	//readmc.DisableAutoScheme();
511	if (!set.AddFilesOn(readmc))
512	return kFALSE;
513
514	readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
515	readmc.EnableBranch("MMcEvtBasic.fEnergy");
516
517	mh1.SetLogy();
518	mh1.SetLogz();
519	mh1.SetName("ThetaE");
520
521	MFillH fill0(&mh1);
522	//fill0.SetDrawOption("projx only");
523
524	MBinning bins2 = (MBinning)plist.FindObject("BinningEnergyEst");
525	MBinning bins3 = (MBinning)plist.FindObject("BinningTheta");
526	if (bins2 && bins3)
527	{
528	bins2->SetName("BinningThetaEY");
529	bins3->SetName("BinningThetaEX");
530	}
531	tlist1.AddToList(&readmc);
532	tlist1.AddToList(&weight);
533
534	temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
535	MH3 mh3mc(temp1);
536
537	MFEventSelector2 sel1(mh3mc);
538	sel1.SetHistIsProbability();
539
540	fill0.SetFilter(&sel1);
541
542	if (!fRawMc)
543	tlist1.AddToList(&sel1);
544	tlist1.AddToList(&fill0);
545
546	// by setting it here it is distributed to all consecutive tasks
547	tlist1.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
548
549	MEvtLoop loop1("IntermediateLoop"); // *** fName ***
550	loop1.SetParList(&plist);
551	loop1.SetLogStream(fLog);
552	loop1.SetDisplay(fDisplay);
553
554	if (!SetupEnv(loop1))
555	return kFALSE;
556
557	if (!loop1.Eventloop(fMaxEvents))
558	{
559	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
560	return kFALSE;
561	}
562
563	if (!loop1.GetDisplay())
564	{
565	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
566	return kFALSE;
567	}
568
569	if (bins2 && bins3)
570	{
571	bins2->SetName("BinningEnergyEst");
572	bins3->SetName("BinningTheta");
573	}
574
575	return kTRUE;
576	}
577
578	TString MJSpectrum::FormString(const TF1 &f, Byte_t type)
579	{
580	const Double_t p0 = f.GetParameter(0);
581	const Double_t p1 = f.GetParameter(1);
582
583	const Double_t e0 = f.GetParError(0);
584	const Double_t e1 = f.GetParError(1);
585
586	const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
587	const Double_t exp = TMath::Power(10, np);
588
589	const Float_t fe0 = TMath::Log10(e0);
590	const Float_t fe1 = TMath::Log10(e1);
591
592	// calc number of (gueltige ziffern) digits to be displayed
593	const char *f0 = fe0-TMath::Floor(fe0)>0.3 ? "3.1" : "1.0";
594	const char *f1 = fe1-TMath::Floor(fe1)>0.3 ? "3.1" : "1.0";
595
596	TString str;
597	switch (type)
598	{
599	case 0:
600	{
601	const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}", f1, f1);
602	const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
603
604	str = Form(fmt0.Data(), p1/exp, e1/exp, np);
605	str += Form(fmt1.Data(), p0, e0);
606	str += "\\frac{ph}{TeVm^{2}s}";
607	}
608	break;
609	case 1:
610	str = Form("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF());
611	break;
612	case 2:
613	str = Form("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(), f.GetNDF()));
614	break;
615	}
616	return str;
617	}
618
619	TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
620	{
621	TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
622	f.SetParameter(0, -2.5);
623	f.SetParameter(1, spectrum.GetBinContent(spectrum.GetXaxis()->FindFixBin(1000)));
624	f.SetLineColor(kBlue);
625	spectrum.Fit(&f, "NI", "", 80, 1150); // M skipped
626	f.DrawCopy("same");
627
628	TString str = FormString(f);
629
630	TLatex tex;
631	tex.SetTextSize(0.045);
632	tex.SetBit(TLatex::kTextNDC);
633	tex.SetTextAlign(31);
634	tex.DrawLatex(0.89, 0.935, str);
635
636	str = FormString(f, 1);
637	tex.DrawLatex(0.89, 0.83, str);
638
639	str = FormString(f, 2);
640	tex.DrawLatex(0.89, 0.735, str);
641
642	TArrayD res(2);
643	res[0] = f.GetParameter(0);
644	res[1] = f.GetParameter(1);
645	return res;
646	}
647
648	// --------------------------------------------------------------------------
649	//
650	// Calculate the final spectrum from:
651	// - collection area
652	// - excess
653	// - correction coefficients
654	// - ontime
655	// and display it
656	//
657	TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
658	{
659	// Create copies of the histograms
660	TH1D collarea(area.GetHEnergy());
661	TH1D spectrum(excess);
662	TH1D weights;
663
664	// Get spill-over corrections from energy estimation
665	hest.GetWeights(weights);
666
667	// Print effective on-time
668	cout << "Effective On time: " << ontime << "s" << endl;
669
670	// Setup spectrum plot
671	spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
672	spectrum.SetYTitle("N/sm^{2}");
673	spectrum.SetDirectory(NULL);
674	spectrum.SetBit(kCanDelete);
675
676	// Divide by collection are and on-time
677	spectrum.Scale(1./ontime);
678	spectrum.Divide(&collarea);
679
680	// Draw spectrum before applying spill-over corrections
681	TCanvas &c1 = fDisplay->AddTab("Spectrum");
682	c1.Divide(2,2);
683	c1.cd(1);
684	gPad->SetBorderMode(0);
685	gPad->SetLogx();
686	gPad->SetLogy();
687	gPad->SetGridx();
688	gPad->SetGridy();
689	collarea.DrawCopy();
690
691	c1.cd(2);
692	gPad->SetBorderMode(0);
693	gPad->SetLogx();
694	gPad->SetLogy();
695	gPad->SetGridx();
696	gPad->SetGridy();
697	TH1D spec=(TH1D)spectrum.DrawCopy();
698	FitSpectrum(*spec);
699
700	c1.cd(3);
701	gPad->SetBorderMode(0);
702	gPad->SetLogx();
703	gPad->SetLogy();
704	gPad->SetGridx();
705	gPad->SetGridy();
706	weights.DrawCopy();
707
708	// Apply spill-over correction (done't take the errors into account)
709	// They are supposed to be identical with the errors of the
710	// collection area and cancel out.
711	//spectrum.Multiply(&weights);
712	spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
713	spectrum.SetBit(TH1::kNoStats);
714
715	for (int i=0; i<excess.GetNbinsX(); i++)
716	{
717	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
718	spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1));
719
720	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
721	spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
722	}
723
724	c1.cd(4);
725	gPad->SetBorderMode(0);
726	gPad->SetLogx();
727	gPad->SetLogy();
728	gPad->SetGridx();
729	gPad->SetGridy();
730	spectrum.SetMinimum(1e-12);
731	spectrum.SetXTitle("E [GeV]");
732	spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
733	spec = (TH1D*)spectrum.DrawCopy();
734
735	// Calculate Spectrum * E^2
736	for (int i=0; i<spectrum.GetNbinsX(); i++)
737	{
738	const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)spectrum.GetBinLowEdge(i+2))1e-3;
739
740	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)ee);
741	spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) ee);
742	}
743
744	// Display dN/dE*E^2 for conveinience
745	fDisplay->AddTab("Check");
746	gPad->SetLogx();
747	gPad->SetLogy();
748	gPad->SetGrid();
749
750	spectrum.SetName("FluxStd");
751	spectrum.SetMarkerStyle(kFullDotMedium);
752	spectrum.SetTitle("Differential flux times E^{2}");
753	spectrum.SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
754	spectrum.SetDirectory(0);
755	spectrum.DrawCopy();
756
757	// Fit Spectrum
758	c1.cd(4);
759	return FitSpectrum(*spec);
760
761	/*
762	TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
763	f.SetParameter(0, -2.87);
764	f.SetParameter(1, 1.9e-6);
765	f.SetLineColor(kGreen);
766	spectrum.Fit(&f, "NIM", "", 100, 5000);
767	f.DrawCopy("same");
768
769	const Double_t p0 = f.GetParameter(0);
770	const Double_t p1 = f.GetParameter(1);
771
772	const Double_t e0 = f.GetParError(0);
773	const Double_t e1 = f.GetParError(1);
774
775	const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
776	const Double_t exp = TMath::Power(10, np);
777
778	TString str;
779	str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
780	str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
781	str += "\\frac{ph}{TeVm^{2}s}";
782
783	TLatex tex;
784	tex.SetTextSize(0.045);
785	tex.SetBit(TLatex::kTextNDC);
786	tex.DrawLatex(0.45, 0.935, str);
787
788	str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
789	tex.DrawLatex(0.70, 0.83, str);
790
791	TArrayD res(2);
792	res[0] = f.GetParameter(0);
793	res[1] = f.GetParameter(1);
794
795	return res;
796	*/
797	/*
798	// Plot other spectra from Whipple
799	f.SetParameter(0, -2.45);
800	f.SetParameter(1, 3.3e-7);
801	f.SetRange(300, 8000);
802	f.SetLineColor(kBlack);
803	f.SetLineStyle(kDashed);
804	f.DrawCopy("same");
805
806	// Plot other spectra from Cangaroo
807	f.SetParameter(0, -2.53);
808	f.SetParameter(1, 2.0e-7);
809	f.SetRange(7000, 50000);
810	f.SetLineColor(kBlack);
811	f.SetLineStyle(kDashed);
812	f.DrawCopy("same");
813
814	// Plot other spectra from Robert
815	f.SetParameter(0, -2.59);
816	f.SetParameter(1, 2.58e-7);
817	f.SetRange(150, 1500);
818	f.SetLineColor(kBlack);
819	f.SetLineStyle(kDashed);
820	f.DrawCopy("same");
821
822	// Plot other spectra from HEGRA
823	f.SetParameter(0, -2.61);
824	f.SetParameter(1, 2.7e-7);
825	f.SetRange(1000, 20000);
826	f.SetLineColor(kBlack);
827	f.SetLineStyle(kDashed);
828	f.DrawCopy("same");
829	*/
830	}
831
832	// --------------------------------------------------------------------------
833	//
834	// Scale some image parameter plots using the scale factor and plot them
835	// together with the corresponding MC histograms.
836	// Called from DisplaySize
837	//
838	Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char name, const char tab, const char *plot, Double_t scale) const
839	{
840	TString same(name);
841	same += "Same";
842
843	TH1 h1 = (TH1)arr.FindObjectInCanvas(name, "TH1F", tab);
844	TH1 h2 = (TH1)arr.FindObjectInCanvas(same, "TH1F", tab);
845	if (!h1 \|\| !h2)
846	return kFALSE;
847
848	TObject *obj = plist.FindObject(plot);
849	if (!obj)
850	{
851	*fLog << warn << plot << " not in parameter list... skipping." << endl;
852	return kFALSE;
853	}
854
855	TH1 h3 = (TH1)obj->FindObject(name);
856	if (!h3)
857	{
858	*fLog << warn << name << " not found in " << plot << "... skipping." << endl;
859	return kFALSE;
860	}
861
862
863	const MAlphaFitter fit = (MAlphaFitter)plist.FindObject("MAlphaFitter");
864	const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
865
866	gPad->SetBorderMode(0);
867	h2->SetLineColor(kBlack);
868	h3->SetLineColor(kBlue);
869	h2->Add(h1, -ascale);
870
871	//h2->Scale(1./ontime); //h2->Integral());
872	h3->Scale(scale); //h3->Integral());
873
874	h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
875
876	h2 = h2->DrawCopy();
877	h3 = h3->DrawCopy("same");
878
879	// Don't do this on the original object!
880	h2->SetStats(kFALSE);
881	h3->SetStats(kFALSE);
882
883	return kTRUE;
884	}
885
886	// --------------------------------------------------------------------------
887	//
888	// Take a lot of histograms and plot them together in one plot.
889	// Calls PlotSame
890	//
891	Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
892	{
893	*fLog << inf << "Reading from file: " << fPathIn << endl;
894
895	TFile file(fPathIn, "READ");
896	if (!file.IsOpen())
897	{
898	*fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
899	return kFALSE;
900	}
901
902	file.cd();
903	MStatusArray arr;
904	if (arr.Read()<=0)
905	{
906	*fLog << "MStatusDisplay not found in file... abort." << endl;
907	return kFALSE;
908	}
909
910	TH1 excess = (TH1D)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
911	if (!excess)
912	return kFALSE;
913
914	// ------------------- Plot excess versus size -------------------
915
916	TCanvas &c = fDisplay->AddTab("Excess");
917	c.Divide(3,2);
918	c.cd(1);
919	gPad->SetBorderMode(0);
920	gPad->SetLogx();
921	gPad->SetLogy();
922	gPad->SetGridx();
923	gPad->SetGridy();
924
925	excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
926	excess = excess->DrawCopy("E2");
927	// Don't do this on the original object!
928	excess->SetStats(kFALSE);
929	excess->SetMarkerStyle(kFullDotMedium);
930	excess->SetFillColor(kBlack);
931	excess->SetFillStyle(0);
932	excess->SetName("Excess ");
933	excess->SetDirectory(0);
934
935	TObject *o=0;
936	if ((o=plist.FindObject("ExcessMC")))
937	{
938	TH1 histsel = (TH1F)o->FindObject("");
939	if (histsel)
940	{
941	if (scale<0)
942	scale = excess->Integral()/histsel->Integral();
943
944	histsel->Scale(scale);
945	histsel->SetLineColor(kBlue);
946	histsel->SetBit(kCanDelete);
947	histsel = histsel->DrawCopy("E1 same");
948	// Don't do this on the original object!
949	histsel->SetStats(kFALSE);
950
951	fLog->Separator("Kolmogorov Test");
952	histsel->KolmogorovTest(excess, "DX");
953	fLog->Separator("Chi^2 Test");
954	const Double_t p = histsel->Chi2Test(excess, "P");
955
956	TLatex tex;
957	tex.SetBit(TLatex::kTextNDC);
958	tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
959	}
960	}
961
962	// -------------- Comparison of Image Parameters --------------
963	c.cd(2);
964	PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale);
965
966	c.cd(3);
967	PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
968
969	c.cd(4);
970	PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale);
971
972	c.cd(5);
973	PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale);
974
975	c.cd(6);
976	PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale);
977
978	return kTRUE;
979	}
980
981	// --------------------------------------------------------------------------
982	//
983	void MJSpectrum::DisplayCutEfficiency(const MHCollectionArea &area0, const MHCollectionArea &area1) const
984	{
985	if (!fDisplay)
986	return;
987
988	const TH1D &trig = area0.GetHEnergy();
989	TH1D &cut = (TH1D&)*area1.GetHEnergy().Clone();
990
991	fDisplay->AddTab("CutEff");
992
993	gPad->SetBorderMode(0);
994	gPad->SetFrameBorderMode(0);
995	gPad->SetLogx();
996	gPad->SetGridx();
997	gPad->SetGridy();
998
999	cut.Divide(&trig);
1000	cut.Scale(100);
1001	cut.SetNameTitle("CutEff", "Background Supression: Cut efficiency (after star)");
1002	cut.SetYTitle("\\eta [%]");
1003	cut.SetDirectory(0);
1004	cut.SetMinimum(0);
1005	cut.SetMaximum(100);
1006	cut.SetBit(kCanDelete);
1007	cut.Draw();
1008
1009	TLine line;
1010	line.SetLineColor(kBlue);
1011	line.SetLineWidth(2);
1012	line.SetLineStyle(kDashed);
1013	line.DrawLine(cut.GetBinLowEdge(1), 50, cut.GetBinLowEdge(cut.GetNbinsX()+1), 50);
1014	}
1015
1016	Bool_t MJSpectrum::Process(const MDataSet &set)
1017	{
1018	if (!set.IsValid())
1019	{
1020	*fLog << err << "ERROR - DataSet invalid!" << endl;
1021	return kFALSE;
1022	}
1023
1024	CheckEnv();
1025
1026	// --------------------------------------------------------------------------------
1027
1028	*fLog << inf;
1029	fLog->Separator(GetDescriptor());
1030	*fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
1031	*fLog << endl;
1032
1033	MBinning bins1("BinningAlpha");
1034	MBinning bins2("BinningEnergyEst");
1035	MBinning bins3("BinningTheta");
1036	MBinning bins4("BinningFalseSource");
1037	MBinning bins5("BinningWidth");
1038	MBinning bins6("BinningLength");
1039	MBinning bins7("BinningDist");
1040	MBinning bins8("BinningMaxDist");
1041	MBinning bins9("BinningM3Long");
1042	MBinning bins0("BinningConc1");
1043
1044	MAlphaFitter fit;
1045
1046	MParList plist;
1047	plist.AddToList(&bins1);
1048	plist.AddToList(&bins3);
1049	plist.AddToList(&bins4);
1050	plist.AddToList(&bins5);
1051	plist.AddToList(&bins6);
1052	plist.AddToList(&bins7);
1053	plist.AddToList(&bins8);
1054	plist.AddToList(&bins9);
1055	plist.AddToList(&bins0);
1056	plist.AddToList(&fit);
1057
1058	TH1D temp1, size;
1059	const Float_t ontime = ReadInput(plist, temp1, size);
1060	if (ontime<0)
1061	{
1062	*fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
1063	return kFALSE;
1064	}
1065
1066	plist.AddToList(&bins2);
1067
1068	// Initialize weighting to a new spectru
1069	// m as defined in the resource file
1070	MMcSpectrumWeight weight;
1071	if (!InitWeighting(set, weight))
1072	return kFALSE;
1073
1074	// Print the setup of the MAlphaFitter
1075	PrintSetup(fit);
1076	bins3.SetEdges(temp1, 'x');
1077
1078	// Read the MC distribution as produced by corsika into
1079	// temp2 (1D) and apply the weights previously determined
1080	TH1D temp2(temp1);
1081	if (!ReadOrigMCDistribution(set, temp2, weight))
1082	return kFALSE;
1083
1084	// Calculate the weights according to the zenith angle distribution
1085	// of the raw-data which have to be applied to the MC events
1086	if (!GetThetaDistribution(temp1, temp2))
1087	return kFALSE;
1088
1089	// Tell the weighting task about the ZA-weights
1090	if (!fNoThetaWeights)
1091	weight.SetWeightsZd(&temp1);
1092
1093	// Refill excess histogram to determin the excess events
1094	TH1D excess;
1095	if (!Refill(plist, excess))
1096	return kFALSE;
1097
1098	TH2D hist;
1099	MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
1100	if (fSimpleMode)
1101	{
1102	// Now we read the MC distribution as produced by corsika again
1103	// with the spectral weights applied into a histogram vs.
1104	// zenith angle and energy
1105	hist.UseCurrentStyle();
1106	MH::SetBinning(&hist, &bins3/temp1.GetXaxis()/, &bins2/excess.GetXaxis()/);
1107	if (!ReadOrigMCDistribution(set, hist, weight))
1108	return kFALSE;
1109
1110	// No we apply the the zenith-angle-weights to the corsika produced
1111	// MC distribution. Unfortunately this must be done manually
1112	// because we are multiplying column by column
1113	if (!fRawMc)
1114	{
1115	for (int y=0; y<hist.GetNbinsY(); y++)
1116	for (int x=0; x<hist.GetNbinsX(); x++)
1117	{
1118	hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
1119	hist.SetBinError(x, y, hist.GetBinError(x, y) *temp1.GetBinContent(x));
1120	}
1121	}
1122	}
1123	else
1124	{
1125	// This rereads the original MC spectrumand aaplies both
1126	// weights spectral weights and ZA-weights.
1127	weight.SetNameMcEvt("MMcEvtBasic");
1128	if (!IntermediateLoop(plist, mh1, temp1, set, weight))
1129	return kFALSE;
1130	weight.SetNameMcEvt();
1131	}
1132
1133	DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1134
1135	// ------------------------- Final loop --------------------------
1136
1137	*fLog << endl;
1138	fLog->Separator("Calculate efficiencies");
1139	*fLog << endl;
1140
1141	MTaskList tlist2;
1142	plist.Replace(&tlist2);
1143
1144	MReadMarsFile read("Events");
1145	read.DisableAutoScheme();
1146	if (!set.AddFilesOn(read))
1147	return kFALSE;
1148
1149	// Selector to get correct (final) theta-distribution
1150	temp1.SetXTitle("MPointingPos.fZd");
1151
1152	MH3 mh3(temp1);
1153
1154	MFEventSelector2 sel2(mh3);
1155	sel2.SetHistIsProbability();
1156
1157	MContinue contsel(&sel2);
1158	contsel.SetInverted();
1159
1160	// Get correct source position
1161	//MSrcPosCalc calc;
1162
1163	// Calculate corresponding Hillas parameters
1164	/*
1165	MHillasCalc hcalc1;
1166	MHillasCalc hcalc2("MHillasCalcAnti");
1167	hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
1168	hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
1169	hcalc2.SetNameHillasSrc("MHillasSrcAnti");
1170	hcalc2.SetNameSrcPosCam("MSrcPosAnti");
1171	*/
1172
1173	// Fill collection area and energy estimator (unfolding)
1174	// Make sure to use the same binning for MHCollectionArea and MHEnergyEst
1175	MHCollectionArea area0("TriggerArea");
1176	MHCollectionArea area1;
1177	area0.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1178	area1.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1179
1180	MHEnergyEst hest;
1181
1182	MFillH fill30(&area0, "", "FillTriggerArea");
1183	MFillH fill31(&area1, "", "FillCollectionArea");
1184	MFillH fill4(&hest, "", "FillEnergyEst");
1185	MFillH fill5("MHThreshold", "", "FillThreshold");
1186	fill30.SetWeight();
1187	fill31.SetWeight();
1188	fill4.SetWeight();
1189	fill5.SetWeight();
1190	fill30.SetNameTab("TrigArea");
1191	fill31.SetNameTab("ColArea");
1192	fill4.SetNameTab("E-Est");
1193	fill5.SetNameTab("Threshold");
1194
1195	MH3 hsize("MHillas.fSize");
1196	hsize.SetName("ExcessMC");
1197	hsize.Sumw2();
1198
1199	MBinning bins(size, "BinningExcessMC");
1200	plist.AddToList(&hsize);
1201	plist.AddToList(&bins);
1202
1203	MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
1204	MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
1205	MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
1206	MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
1207	MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
1208	MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
1209	MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
1210	MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC");
1211	fill1a.SetNameTab("PreCut");
1212	fill2a.SetNameTab("PostCut");
1213	fill3a.SetNameTab("VsSize");
1214	fill4a.SetNameTab("HilExt");
1215	fill5a.SetNameTab("HilSrc");
1216	fill6a.SetNameTab("ImgPar");
1217	fill7a.SetNameTab("NewPar");
1218	fill8a.SetBit(MFillH::kDoNotDisplay);
1219	fill1a.SetWeight();
1220	fill2a.SetWeight();
1221	fill3a.SetWeight();
1222	fill4a.SetWeight();
1223	fill5a.SetWeight();
1224	fill6a.SetWeight();
1225	fill7a.SetWeight();
1226	fill8a.SetWeight();
1227
1228	MTaskEnv taskenv0("CalcHadronness");
1229	taskenv0.SetDefault(fCalcHadronness);
1230
1231	MEnergyEstimate est;
1232	MTaskEnv taskenv1("EstimateEnergy");
1233	taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
1234
1235	tlist2.AddToList(&read);
1236	// If no weighting should be applied but the events should
1237	// be thrown away according to the theta distribution
1238	// it is enabled here
1239	if (!fRawMc && fNoThetaWeights)
1240	tlist2.AddToList(&contsel);
1241	//tlist2.AddToList(&calc);
1242	//tlist2.AddToList(&hcalc1);
1243	//tlist2.AddToList(&hcalc2);
1244	tlist2.AddToList(&weight);
1245	tlist2.AddToList(&fill1a);
1246	tlist2.AddToList(&fill30);
1247	tlist2.AddToList(fCut0);
1248	tlist2.AddToList(&taskenv0);
1249	tlist2.AddToList(fCut1);
1250	tlist2.AddToList(fCut2);
1251	tlist2.AddToList(fCut3);
1252	tlist2.AddToList(&taskenv1);
1253	tlist2.AddToList(&fill31);
1254	tlist2.AddToList(&fill4);
1255	tlist2.AddToList(&fill5);
1256	tlist2.AddToList(&fill2a);
1257	tlist2.AddToList(&fill3a);
1258	tlist2.AddToList(&fill4a);
1259	tlist2.AddToList(&fill5a);
1260	tlist2.AddToList(&fill6a);
1261	tlist2.AddToList(&fill7a);
1262	tlist2.AddToList(&fill8a);
1263	//tlist2.AddToList(&fill9a);
1264
1265	// by setting it here it is distributed to all consecutive tasks
1266	tlist2.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
1267
1268	MEvtLoop loop2("FillMonteCarlo"); // *** fName ***
1269	loop2.SetParList(&plist);
1270	loop2.SetDisplay(fDisplay);
1271	loop2.SetLogStream(fLog);
1272
1273	if (!SetupEnv(loop2))
1274	return kFALSE;
1275
1276	if (!loop2.Eventloop(fMaxEvents))
1277	{
1278	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
1279	return kFALSE;
1280	}
1281
1282	if (!loop2.GetDisplay())
1283	{
1284	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
1285	return kFALSE;
1286	}
1287
1288	gLog.Separator("Energy Estimator");
1289	if (plist.FindObject("EstimateEnergy"))
1290	plist.FindObject("EstimateEnergy")->Print();
1291
1292	gLog.Separator("Spectrum");
1293
1294	// -------------------------- Spectrum ----------------------------
1295
1296	// Calculate and display spectrum (N/TeVsm^2 at 1TeV)
1297	TArrayD res(DisplaySpectrum(area1, excess, hest, ontime));
1298
1299	// Spectrum fitted (convert res[1] from TeV to GeV)
1300	TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
1301
1302	// Number of events this spectrum would produce per s and m^2
1303	Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
1304
1305	// scale with effective collection area to get the event rate (N/s)
1306	// scale with the effective observation time to absolute observed events
1307	n = area1.GetCollectionAreaAbs()ontime; // N
1308
1309	// Now calculate the scale factor from the number of events
1310	// produced and the number of events which should have been
1311	// observed with our telescope in the time ontime
1312	const Double_t scale = n/area1.GetEntries();
1313
1314	// Print normalization constant
1315	cout << "MC normalization factor: " << scale << endl;
1316
1317	// Display cut efficiency
1318	DisplayCutEfficiency(area0, area1);
1319
1320	// Overlay normalized plots
1321	DisplaySize(plist, scale);
1322
1323	// check if output should be written
1324	if (fPathOut.IsNull())
1325	return kTRUE;
1326
1327	// Write the output
1328	TObjArray cont;
1329	cont.Add((TObject*)GetEnv());
1330
1331	if (fDisplay)
1332	cont.Add(fDisplay);
1333
1334	return WriteContainer(cont, "", "RECREATE");
1335	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: