Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6957) +++ trunk/MagicSoft/Mars/Changelog (revision 6958) @@ -29,4 +29,34 @@ * mimage/MHVsSize.[h,cc]: - added more displayed variables + + * mjobs/MJSpectrum.[h,cc]: + - added + + * mjobs/Makefile, mjobs/JobsLinkDef.h: + - added MJSpectrum + + * Makefile: + - added sponde.cc + + * sponde.cc: + - added + + * mhbase/MBinning.h: + - added Print() to //*MENU* + + * mhflux/MAlphaFitter.cc: + - do not set inf-flag in Print() + + * mhflux/MHAlpha.[h,cc]: + - set all-flag for printing MAlphaFitter("result") + - added Getter for excess histogram + + * mjobs/MDataSet.[h,cc]: + - allow to initialize a TChain from this dataset + + * mjobs/MJCut.cc: + - moved energy estimator back before Cut0 otherwise the + EnergyEst container is missing in summary file + Index: trunk/MagicSoft/Mars/Makefile =================================================================== --- trunk/MagicSoft/Mars/Makefile (revision 6957) +++ trunk/MagicSoft/Mars/Makefile (revision 6958) @@ -20,5 +20,5 @@ # #PROGRAMS = readraw merpp mars test mona status -PROGRAMS = readdaq merpp readraw sinope callisto star ganymed showlog showplot mars +PROGRAMS = readdaq merpp readraw sinope callisto star ganymed sponde showlog showplot mars SOLIB = libmars.so CINT = M Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 6957) +++ trunk/MagicSoft/Mars/NEWS (revision 6958) @@ -1,4 +1,7 @@ -*-*- END -*-*- *** Version + + - added a full featured spectrum program (sponde) which reads + ganymed output and monte carlos and compiles a spectrum Index: trunk/MagicSoft/Mars/mhbase/MBinning.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MBinning.h (revision 6957) +++ trunk/MagicSoft/Mars/mhbase/MBinning.h (revision 6958) @@ -108,5 +108,5 @@ void Apply(TH1 &) const; - void Print(Option_t *o="") const; + void Print(Option_t *o="") const; //*MENU* Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 6957) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 6958) @@ -328,5 +328,5 @@ void MAlphaFitter::Print(Option_t *o) const { - *fLog << inf << GetDescriptor() << ": Fitting..." << endl; + *fLog << GetDescriptor() << ": Fitting..." << endl; *fLog << " ...background from " << fBgMin << " to " << fBgMax << endl; *fLog << " ...signal to " << fSigMax << " (integrate into bin at " << fSigInt << ")" << endl; Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 6957) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 6958) @@ -352,4 +352,5 @@ fit->Copy(fFit); + *fLog << inf; fFit.Print(); @@ -864,4 +865,5 @@ // Store the final result in fFit + *fLog << all; fFit.Print("result"); Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 6957) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 6958) @@ -137,4 +137,6 @@ const MAlphaFitter &GetAlphaFitter() const { return fFit; } + const TH1D &GetHEnergy() const { return fHEnergy; } + void SetOffData(const MHAlpha &h) { Index: trunk/MagicSoft/Mars/mjobs/JobsLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/JobsLinkDef.h (revision 6957) +++ trunk/MagicSoft/Mars/mjobs/JobsLinkDef.h (revision 6958) @@ -18,4 +18,5 @@ #pragma link C++ class MJCut+; #pragma link C++ class MJOptimize+; +#pragma link C++ class MJSpectrum+; #endif Index: trunk/MagicSoft/Mars/mjobs/MDataSet.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MDataSet.cc (revision 6957) +++ trunk/MagicSoft/Mars/mjobs/MDataSet.cc (revision 6958) @@ -71,4 +71,5 @@ #include +#include #include #include // TSystem::ExpandPath @@ -234,8 +235,6 @@ } -Bool_t MDataSet::AddSequencesToList(const TList &list, MRead &read) -{ - MDirIter files; - +Bool_t MDataSet::AddSequencesFromList(const TList &list, MDirIter &files) +{ TIter Next(const_cast(&list)); TObject *o=0; @@ -258,5 +257,20 @@ files.Print(); } - + return kTRUE; +} + +Bool_t MDataSet::AddFilesOn(MRead &read) const +{ + MDirIter files; + if (!AddSequencesFromList(fSequencesOn, files)) + return kFALSE; + return read.AddFiles(files)>0; +} + +Bool_t MDataSet::AddFilesOff(MRead &read) const +{ + MDirIter files; + if (!AddSequencesFromList(fSequencesOff, files)) + return kFALSE; return read.AddFiles(files)>0; } @@ -269,12 +283,42 @@ } -Bool_t MDataSet::AddFilesOn(MRead &read) const -{ - return AddSequencesToList(fSequencesOn, read); -} - -Bool_t MDataSet::AddFilesOff(MRead &read) const -{ - return AddSequencesToList(fSequencesOff, read); +Int_t MDataSet::AddFilesToChain(MDirIter &files, TChain &chain) +{ + Int_t num=0; + while (1) + { + const TString fname = files.Next(); + if (fname.IsNull()) + break; + + const Int_t n = chain.Add(fname); + if (n<=0) + return kFALSE; + num += n; + } + return num; +} + +Bool_t MDataSet::AddFilesOn(TChain &chain) const +{ + MDirIter files; + if (!AddSequencesFromList(fSequencesOn, files)) + return kFALSE; + return AddFilesToChain(files, chain)>0; +} + +Bool_t MDataSet::AddFilesOff(TChain &chain) const +{ + MDirIter files; + if (!AddSequencesFromList(fSequencesOff, files)) + return kFALSE; + return AddFilesToChain(files, chain)>0; +} + +Bool_t MDataSet::AddFiles(TChain &read) const +{ + const Bool_t rc1 = AddFilesOff(read); + const Bool_t rc2 = AddFilesOn(read); + return rc1 && rc2; } Index: trunk/MagicSoft/Mars/mjobs/MDataSet.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MDataSet.h (revision 6957) +++ trunk/MagicSoft/Mars/mjobs/MDataSet.h (revision 6958) @@ -10,5 +10,8 @@ #endif +class TChain; + class MRead; +class MDirIter; class MPointingPos; @@ -54,9 +57,14 @@ // Setter - static Bool_t AddSequencesToList(const TList &list, MRead &read); + static Bool_t AddSequencesFromList(const TList &list, MDirIter &files); + static Int_t AddFilesToChain(MDirIter &files, TChain &chain); Bool_t AddFiles(MRead &read) const; Bool_t AddFilesOn(MRead &read) const; Bool_t AddFilesOff(MRead &read) const; + + Bool_t AddFiles(TChain &read) const; + Bool_t AddFilesOn(TChain &read) const; + Bool_t AddFilesOff(TChain &read) const; // TObject Index: trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc (revision 6957) +++ trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc (revision 6958) @@ -373,5 +373,5 @@ interlacedcont.Add(&hpndiod); if (fIsRelTimesUpdate) - interlacedcont.Add(&hrelcam); + interlacedcont.Add(&hrelcam); } @@ -601,11 +601,10 @@ // Intensity Calibration histogramming // - MFillH filpin("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode", "FillPINDiode"); - MFillH filbnd("MHCalibrationChargeBlindCam", "MExtractedSignalBlindPixel", "FillBlindCam"); - MFillH filcam("MHCalibrationChargeCam", "MExtractedSignalCam", "FillChargeCam"); - MFillH filtme("MHCalibrationRelTimeCam", "MArrivalTimeCam", "FillRelTime"); + MFillH filpin(&hpndiod, "MExtractedSignalPINDiode", "FillPINDiode"); + MFillH filbnd(&hbndcam, "MExtractedSignalBlindPixel", "FillBlindCam"); + MFillH filcam(&hchacam, "MExtractedSignalCam", "FillChargeCam"); + MFillH filtme(&hrelcam, "MArrivalTimeCam", "FillRelTime"); MFillH filhil("MHCalibrationHiLoCam", "MExtractedSignalCam", "FillHiLoRatio"); MFillH filpul("MHCalibrationPulseTimeCam", "MRawEvtData", "FillPulseTime"); - filpin.SetBit(MFillH::kDoNotDisplay); filbnd.SetBit(MFillH::kDoNotDisplay); Index: trunk/MagicSoft/Mars/mjobs/MJCut.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 6957) +++ trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 6958) @@ -495,4 +495,5 @@ if (fIsWobble) tlist2.AddToList(&hcalc2); + tlist2.AddToList(&taskenv1); tlist2.AddToList(&taskenv2); tlist2.AddToList(&cont0); @@ -502,5 +503,4 @@ tlist2.AddToList(&fill1a); tlist2.AddToList(&cont1); - tlist2.AddToList(&taskenv1); if (!fWriteOnly && !fIsWobble) tlist2.AddToList(&ffs); Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 6958) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 6958) @@ -0,0 +1,726 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 4/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MJSpectrum +// +// Program to calculate spectrum +// +///////////////////////////////////////////////////////////////////////////// +#include "MJSpectrum.h" + +// Root +#include +#include +#include +#include +#include +#include +#include + +// Environment +#include "MLog.h" +#include "MLogManip.h" + +#include "MStatusArray.h" +#include "MStatusDisplay.h" + +// Container +#include "MH3.h" +#include "MBinning.h" +#include "MDataSet.h" + +// Spectrum +#include "../mhflux/MAlphaFitter.h" +#include "../mhflux/MHAlpha.h" +#include "../mhflux/MHCollectionArea.h" +#include "../mhflux/MHEnergyEst.h" + +// Eventloop +#include "MEvtLoop.h" +#include "MTaskList.h" +#include "MParList.h" + +// Tasks/Filter +#include "MReadMarsFile.h" +#include "MFEventSelector2.h" +#include "MFDataMember.h" +#include "MFillH.h" +#include "MHillasCalc.h" +#include "MSrcPosCalc.h" +#include "MContinue.h" + +ClassImp(MJSpectrum); + +using namespace std; + +MJSpectrum::MJSpectrum(const char *name, const char *title) + : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0), + fRefill(kFALSE), fSimpleMode(kTRUE) +{ + fName = name ? name : "MJSpectrum"; + fTitle = title ? title : "Standard program to calculate spectrum"; +} + +MJSpectrum::~MJSpectrum() +{ + if (fCut0) + delete fCut0; + if (fCut1) + delete fCut1; + if (fCut2) + delete fCut2; + if (fCut3) + delete fCut3; + if (fEstimateEnergy) + delete fEstimateEnergy; +} + +Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const +{ + if (task) + { + delete task; + task = 0; + } + + task = (MTask*)gFile->Get(name); + if (!task) + { + *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl; + return kFALSE; + } + if (!task->InheritsFrom(MTask::Class())) + { + *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl; + delete task; + return kFALSE; + } + + task->SetName(name); + + return kTRUE; +} + +Bool_t MJSpectrum::ReadInput(const MParList &plist) +{ + const TString fname = fPathIn; + + *fLog << inf << "Reading from file: " << fname << endl; + + TFile file(fname, "READ"); + if (!file.IsOpen()) + { + *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl; + return kFALSE; + } + + if (!ReadTask(fCut0, "Cut0")) + return kFALSE; + if (!ReadTask(fCut1, "Cut1")) + return kFALSE; + if (!ReadTask(fCut2, "Cut2")) + return kFALSE; + if (!ReadTask(fCut3, "Cut3")) + return kFALSE; + + if (!ReadTask(fEstimateEnergy, "EstimateEnergy")) + return kFALSE; + + TObjArray arr; + + TIter Next(plist); + TObject *o=0; + while ((o=Next())) + if (o->InheritsFrom(MBinning::Class())) + arr.Add(o); + + arr.Add(plist.FindObject("MAlphaFitter")); + + return ReadContainer(arr); +} + +void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const +{ + gLog.Separator("Alpha Fitter"); + *fLog << all; + fit.Print(); + + gLog.Separator("Used Cuts"); + fCut0->Print(); + fCut1->Print(); + fCut2->Print(); + fCut3->Print(); + + gLog.Separator("Energy Estimator"); + fEstimateEnergy->Print(); +} + +Float_t MJSpectrum::ReadHistograms(TH1D &h1, TH1D &h2) const +{ + TFile file(fPathIn, "READ"); + if (!file.IsOpen()) + { + *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl; + return -1; + } + + MStatusArray *arr = (MStatusArray*)file.Get("MStatusDisplay"); + if (!arr) + { + gLog << "MStatusDisplay not found in file... abort." << endl; + return -1; + } + + TH1D *vstime = (TH1D*)arr->FindObjectInCanvas("Theta", "TH1D", "OnTime"); + TH1D *excen = (TH1D*)arr->FindObjectInCanvas("ExcessEnergy", "TH1D", "MHAlpha"); + if (!vstime || !excen) + return -1; + + vstime->Copy(h1); + excen->Copy(h2); + h1.SetDirectory(0); + h2.SetDirectory(0); + + if (fDisplay) + arr->DisplayIn(*fDisplay, "MHAlpha"); + + delete arr; + + return vstime->Integral(); +} + +Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const +{ + // Some debug output + fLog->Separator("Compiling original MC distribution"); + + *fLog << inf << "Please stand by, this may take a while..." << flush; + + if (fDisplay) + fDisplay->SetStatusLine1("Compiling MC distribution..."); + + // Create chain + TChain chain("OriginalMC"); + set.AddFilesOn(chain); + + // Prepare histogram + h.Reset(); + + // Fill histogram from chain + h.SetDirectory(gROOT); + if (h.InheritsFrom(TH2::Class())) + { + h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)"); + h.SetXTitle("\\Theta [\\circ]"); + h.SetYTitle("E [GeV]"); + h.SetZTitle("Counts"); + chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff"); + } + else + { + h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)"); + h.SetXTitle("\\Theta [\\circ]"); + h.SetYTitle("Counts"); + chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff"); + } + h.SetDirectory(0); + + *fLog << "done." << endl; + if (fDisplay) + fDisplay->SetStatusLine2("done."); + + if (h.GetEntries()>0) + return kTRUE; + + *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl; + + return h.GetEntries()>0; +} + +Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const +{ + // Display some stuff + if (fDisplay) + { + TCanvas &c = fDisplay->AddTab("ZdDist"); + c.Divide(2,2); + + // On-Time vs. Theta + c.cd(1); + gPad->SetBorderMode(0); + temp1.DrawCopy(); + + // Number of MC events (produced) vs Theta + c.cd(2); + gPad->SetBorderMode(0); + temp2.SetName("NVsTheta"); + temp2.DrawCopy(); + + c.cd(4); + gPad->SetBorderMode(0); + + c.cd(3); + gPad->SetBorderMode(0); + } + + // Calculate the Probability + temp1.Divide(&temp2); + temp1.Scale(1./temp1.GetMaximum()); + + // Some cosmetics: Name, Axis, etc. + temp1.SetName("ProbVsTheta"); + temp1.SetTitle("Probability vs. Zenith Angle to choose MC events"); + temp1.SetYTitle("Probability"); + if (fDisplay) + temp1.DrawCopy(); + + return kTRUE; +} + +void MJSpectrum::DisplayResult(const MH3 &mh1) const +{ + if (!fDisplay->CdCanvas("ZdDist")) + return; + + const TH2D &h2 = (const TH2D&)mh1.GetHist(); + + TH1D &proj = *h2.ProjectionX(); + proj.SetNameTitle("ThetaFinal", "Final Theta Distribution"); + proj.SetXTitle("\\Theta [\\circ]"); + proj.SetYTitle("Counts"); + proj.SetLineColor(kBlue); + proj.SetDirectory(0); + proj.SetBit(kCanDelete); + + TVirtualPad *pad = gPad; + + pad->cd(4); + proj.DrawCopy(); + + pad->cd(1); + TH1D *theta = (TH1D*)gPad->FindObject("Theta"); + if (theta) + { + proj.Scale(theta->Integral()/proj.Integral()); + theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum())); + } + proj.Draw("same"); +} + +Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const +{ + *fLog << endl; + fLog->Separator("Refill Excess"); + *fLog << endl; + + MTaskList tlist; + plist.AddToList(&tlist); + + MReadTree read("Events"); + read.DisableAutoScheme(); + read.AddFile(fPathIn); + + MFillH fill1("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha"); + MFillH fill2("MHAlpha", "MHillasSrc", "FillAlpha"); + + MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData"); + MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData"); + + fill1.SetFilter(&f0); + fill2.SetFilter(&f1); + + tlist.AddToList(&read); + tlist.AddToList(fEstimateEnergy); + tlist.AddToList(&f0); + tlist.AddToList(&f1); + tlist.AddToList(&fill1); + tlist.AddToList(&fill2); + + MEvtLoop loop(fName); + loop.SetParList(&plist); + loop.SetDisplay(fDisplay); + loop.SetLogStream(fLog); + + if (!SetupEnv(loop)) + return kFALSE; + + if (!loop.Eventloop()) + { + *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl; + return kFALSE; + } + + tlist.PrintStatistics(); + + if (!loop.GetDisplay()) + { + *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl; + return kFALSE; + } + + const MHAlpha *halpha = (MHAlpha *)plist.FindObject("MHAlpha"); + if (!halpha) + { + *fLog << err << GetDescriptor() << ": MHAlpha not found... abort." << endl; + return kFALSE; + } + + halpha->GetHEnergy().Copy(h2); + h2.SetDirectory(0); + + return kTRUE; +} + +Bool_t MJSpectrum::Process(const MDataSet &set) +{ + if (!set.IsValid()) + { + *fLog << err << "ERROR - DataSet invalid!" << endl; + return kFALSE; + } + + CheckEnv(); + + // -------------------------------------------------------------------------------- + + *fLog << inf; + fLog->Separator(GetDescriptor()); + *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl; + *fLog << endl; + + MBinning bins1("BinningAlpha"); + MBinning bins2("BinningEnergyEst"); + MBinning bins3("BinningTheta"); + MBinning bins4("BinningFalseSource"); + MBinning bins5("BinningWidth"); + MBinning bins6("BinningLength"); + MBinning bins7("BinningDist"); + MBinning bins8("BinningMaxDist"); + + MAlphaFitter fit; + + MParList plist; + plist.AddToList(&bins1); + plist.AddToList(&bins2); + plist.AddToList(&bins3); + plist.AddToList(&bins4); + plist.AddToList(&bins5); + plist.AddToList(&bins6); + plist.AddToList(&bins7); + plist.AddToList(&bins8); + plist.AddToList(&fit); + + if (!ReadInput(plist)) + return kFALSE; + + PrintSetup(fit); + + TH1D temp1, excess; + const Float_t ontime = ReadHistograms(temp1, excess); + if (ontime<0) + { + *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl; + return kFALSE; + } +/* + TH2D hist; + if (fSimpleMode) + { + hist.UseCurrentStyle(); + MH::SetBinning(&hist, temp1.GetXaxis(), excess.GetXaxis()); + if (!ReadOrigMCDistribution(set, hist)) + return kFALSE; + + TCanvas &c =fDisplay->AddTab("test"); + c.Divide(2,2); + c.cd(1); + gPad->SetLogy(); + hist.DrawCopy(); + c.cd(2); + hist.ProjectionX()->Draw(); + c.cd(3); + hist.ProjectionY()->Draw(); + + MH3 mh1(hist); + mh1.SetLogy(); + mh1.SetLogz(); + mh1.SetName("ThetaE"); + } + */ + TH1D temp2(temp1); + if (!ReadOrigMCDistribution(set, temp2)) + return kFALSE; + + if (!GetThetaDistribution(temp1, temp2)) + return kFALSE; + + if (fRefill && !Refill(plist, temp2)) + return kFALSE; + + // ------------------- Loop1 ---------------------------- + + MTaskList tlist1; + plist.Replace(&tlist1); + + MReadMarsFile readmc("OriginalMC"); + //readmc.DisableAutoScheme(); + set.AddFilesOn(readmc); + readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta"); + readmc.EnableBranch("MMcEvtBasic.fEnergy"); + + MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy"); + mh1.SetLogy(); + mh1.SetLogz(); + mh1.SetName("ThetaE"); + + MFillH fill0(&mh1); + //fill0.SetDrawOption("projx only"); + + MBinning binsx(bins3, "BinningThetaEX"); + MBinning binsy(bins2, "BinningThetaEY"); + plist.AddToList(&binsx); + plist.AddToList(&binsy); + tlist1.AddToList(&readmc); + + temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg"); + MH3 mh3mc(temp1); + + MFEventSelector2 sel1(mh3mc); + sel1.SetHistIsProbability(); + + fill0.SetFilter(&sel1); + + tlist1.AddToList(&sel1); + tlist1.AddToList(&fill0); + + MEvtLoop loop1(fName); + loop1.SetParList(&plist); + loop1.SetLogStream(fLog); + loop1.SetDisplay(fDisplay); + + if (!SetupEnv(loop1)) + return kFALSE; + + if (!loop1.Eventloop(fMaxEvents)) + { + *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl; + return kFALSE; + } + + tlist1.PrintStatistics(); + + if (!loop1.GetDisplay()) + { + *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl; + return kFALSE; + } + + DisplayResult(mh1); + + // ------------------------- Final loop -------------------------- + + *fLog << endl; + fLog->Separator("Calculate efficiencies"); + *fLog << endl; + + MTaskList tlist2; + plist.Replace(&tlist2); + + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + set.AddFilesOn(read); + + // Selector to get correct (final) theta-distribution + temp1.SetXTitle("MPointingPos.fZd"); + + MH3 mh3(temp1); + + MFEventSelector2 sel2(mh3); + sel2.SetHistIsProbability(); + + MContinue contsel(&sel2); + contsel.SetInverted(); + + // Get correct source position + MSrcPosCalc calc; + + // Calculate corresponding Hillas parameters + MHillasCalc hcalc1; + MHillasCalc hcalc2("MHillasCalcAnti"); + hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc); + hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc); + hcalc2.SetNameHillasSrc("MHillasSrcAnti"); + hcalc2.SetNameSrcPosCam("MSrcPosAnti"); + + // Fill collection area and energy estimator (unfolding) + // Make sure to use the same binning for MHCollectionArea and MHEnergyEst + MHCollectionArea area; + area.SetHistAll((TH2D&)mh1.GetHist()); + MHEnergyEst hest; + + MFillH fill3(&area); + MFillH fill4(&hest); + + MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre"); + MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost"); + MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost"); + MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost"); + MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost"); + MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost"); + MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost"); + fill1a.SetNameTab("PreCut"); + fill2a.SetNameTab("PostCut"); + fill3a.SetNameTab("VsSize"); + fill4a.SetNameTab("HilExt"); + fill5a.SetNameTab("HilSrc"); + fill6a.SetNameTab("ImgPar"); + fill7a.SetNameTab("NewPar"); + + tlist2.AddToList(&read); + tlist2.AddToList(&contsel); + tlist2.AddToList(&calc); + tlist2.AddToList(&hcalc1); + tlist2.AddToList(&hcalc2); + tlist2.AddToList(&fill1a); + tlist2.AddToList(fCut0); + tlist2.AddToList(fCut1); + tlist2.AddToList(fCut2); + tlist2.AddToList(fCut3); + tlist2.AddToList(fEstimateEnergy); + tlist2.AddToList(&fill3); + tlist2.AddToList(&fill4); + tlist2.AddToList(&fill2a); + tlist2.AddToList(&fill3a); + tlist2.AddToList(&fill4a); + tlist2.AddToList(&fill5a); + tlist2.AddToList(&fill6a); + tlist2.AddToList(&fill7a); + + MEvtLoop loop2(fName); + loop2.SetParList(&plist); + loop2.SetDisplay(fDisplay); + loop2.SetLogStream(fLog); + + if (!SetupEnv(loop2)) + return kFALSE; + + if (!loop2.Eventloop(fMaxEvents)) + { + *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl; + return kFALSE; + } + + tlist2.PrintStatistics(); + + if (!loop2.GetDisplay()) + { + *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl; + return kFALSE; + } + + // -------------------------- Spectrum ---------------------------- + + TH1D collarea(area.GetHEnergy()); + TH1D weights; + hest.GetWeights(weights); + + cout << "Effective On time: " << ontime << "s" << endl; + + excess.SetDirectory(NULL); + excess.SetBit(kCanDelete); + excess.Scale(1./ontime); + excess.Divide(&collarea); + excess.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)"); + excess.SetYTitle("N/sm^{2}"); + + TCanvas &c = fDisplay->AddTab("Spectrum"); + c.Divide(2,2); + c.cd(1); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + collarea.DrawCopy(); + + c.cd(2); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + excess.DrawCopy(); + + c.cd(3); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + weights.DrawCopy(); + + excess.Divide(&weights); + //excess.Multiply(&weights); + excess.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)"); + + for (int i=0; iSetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + excess.DrawCopy(); + + TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4); + f.SetParameter(0, -2.87); + f.SetParameter(1, 1.9e-6); + f.SetLineColor(kGreen); + excess.Fit(&f, "NI", "", 55, 2e4); + f.DrawCopy("same"); + + if (!fPathOut.IsNull()) + fDisplay->SaveAsRoot(fPathOut); + + /* + TString str; + str += "(1.68#pm0.15)10^{-7}"; + str += "(\\frac{E}{TeV})^{-2.59#pm0.06}"; + str += "\\frac{ph}{TeVm^{2}s}"; + + TLatex tex; + tex.DrawLatex(2e2, 7e-5, str); + */ + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 6958) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 6958) @@ -0,0 +1,50 @@ +#ifndef MARS_MJSpectrum +#define MARS_MJSpectrum + +#ifndef MARS_MJob +#include +#endif + +class TH1; +class TH1D; + +class MH3; +class MTask; +class MParList; +class MDataSet; +class MAlphaFitter; + +class MJSpectrum : public MJob +{ +private: + MTask *fCut0; + MTask *fCut1; + MTask *fCut2; + MTask *fCut3; + MTask *fEstimateEnergy; + + Bool_t fRefill; + Bool_t fSimpleMode; + + Bool_t ReadTask(MTask* &task, const char *name) const; + Bool_t ReadInput(const MParList &plist); + void PrintSetup(const MAlphaFitter &fit) const; + Float_t ReadHistograms(TH1D &h1, TH1D &h2) const; + Bool_t ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const; + Bool_t GetThetaDistribution(TH1D &temp1, TH1D &temp2) const; + void DisplayResult(const MH3 &mh1) const; + Bool_t Refill(MParList &plist, TH1D &h) const; + +public: + MJSpectrum(const char *name=NULL, const char *title=NULL); + ~MJSpectrum(); + + Bool_t Process(const MDataSet &set); + + void EnableRefilling(Bool_t b=kTRUE) { fRefill=b; } + void EnableSimpleMode(Bool_t b=kTRUE) { fSimpleMode=b; } + + ClassDef(MJSpectrum, 0) // Proh'gram to calculate spectrum +}; + +#endif Index: trunk/MagicSoft/Mars/mjobs/Makefile =================================================================== --- trunk/MagicSoft/Mars/mjobs/Makefile (revision 6957) +++ trunk/MagicSoft/Mars/mjobs/Makefile (revision 6958) @@ -35,5 +35,6 @@ MJStar.cc \ MJCut.cc \ - MJOptimize.cc + MJOptimize.cc \ + MJSpectrum.cc ############################################################ Index: trunk/MagicSoft/Mars/sponde.cc =================================================================== --- trunk/MagicSoft/Mars/sponde.cc (revision 6958) +++ trunk/MagicSoft/Mars/sponde.cc (revision 6958) @@ -0,0 +1,249 @@ +#include +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MEnv.h" +#include "MArgs.h" +#include "MArray.h" +#include "MDirIter.h" + +#include "MStatusDisplay.h" + +#include "MDataSet.h" +#include "MJSpectrum.h" + +using namespace std; + +static void StartUpMessage() +{ + gLog << all << endl; + + // 1 2 3 4 5 + // 12345678901234567890123456789012345678901234567890123456 + gLog << "========================================================" << endl; + gLog << " Sponde - MARS V" << MARSVER << endl; + gLog << " MARS -- SPectrum ON DEmand" << endl; + gLog << " Compiled on <" << __DATE__ << ">" << endl; + gLog << " Using ROOT v" << ROOTVER << endl; + gLog << "========================================================" << endl; + gLog << endl; +} + +static void Usage() +{ + // 1 2 3 4 5 6 7 8 + // 12345678901234567890123456789012345678901234567890123456789012345678901234567890 + gLog << all << endl; + gLog << "Sorry the usage is:" << endl; + gLog << " sponde [options] inputfile.root mcdataset.txt [outputfile.root]" << endl << endl; + gLog << " Arguments:" << endl; + gLog << " inputfile.root: Ganymed output [and result] file" << endl; + gLog << " mcdataset.txt: Dataset describing the Monte Carlos to be used" << endl; + gLog << " For more details see MDataSet." << endl; + gLog << " outputfile.root: Optional file in which the result is stored" << endl << endl; + gLog << " Root Options:" << endl; + gLog << " -b Batch mode (no graphical output to screen)" << endl<" << endl; + gLog << " --debug-env[=1] Display untouched resources after program execution" << endl; + gLog << " --debug-env=2 Display untouched resources after eventloop setup" << endl; + gLog << " --debug-env=3 Debug setting resources from resource file" << endl; + gLog << " --debug-mem Debug memory usage" << endl << endl; + gLog << endl; + gLog << " -q Quit when job is finished" << endl; + gLog << " -f Force overwrite of existing files" << endl; + gLog << " --print-ds Print Dataset information" << endl; + gLog << " --print-files Print Files taken from Sequences ('+' found, '-' missing)" << endl; + gLog << " --config=sponde.rc Resource file [default=sponde.rc]" << endl; + gLog << endl; + gLog << " --version, -V Show startup message with version number" << endl; + gLog << " -?, -h, --help This help" << endl << endl; + gLog << "Background:" << endl; + gLog << " Sponde is a natural satellite of Jupiter. It was discovered by a team" << endl; + gLog << " of astronomers from the University of Hawaii led by Scott S. Sheppard," << endl; + gLog << " et al in 2001, and given the temporary designation S/2001 J 5. Sponde" << endl; + gLog << " is about 2 kilometers in diameter, and orbits Jupiter at an average" << endl; + gLog << " of 23,487,000 kilometers. It is also designated as Jupiter XXXVI." << endl; + gLog << " It is named after one of the Horae (Hours), which presided over the" << endl; + gLog << " seventh hour (libations poured after lunch). The Hours, godesses of" << endl; + gLog << " the time of day but also of the seasons, were daughters of Zeus and" << endl; + gLog << " Themis." << endl; + gLog << " It belongs to the Pasipha‰ group, irregular retrograde moons orbiting" << endl; + gLog << " Jupiter at distances ranging between 22.8 and 24.1 Gm, and with" << endl; + gLog << " inclinations ranging between 144.5ø and 158.3ø." << endl << endl; +} + +int main(int argc, char **argv) +{ + StartUpMessage(); + + // + // Evaluate arguments + // + MArgs arg(argc, argv, kTRUE); + + if (arg.HasOnly("-V") || arg.HasOnly("--version")) + return 0; + + if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) + { + Usage(); + return 2; + } + + gLog.Setup(arg); + + const TString kConfig = arg.GetStringAndRemove("--config=", "sponde.rc"); + + const Bool_t kPrintSeq = arg.HasOnlyAndRemove("--print-ds"); + const Bool_t kPrintFiles = arg.HasOnlyAndRemove("--print-files"); + const Bool_t kDebugMem = arg.HasOnlyAndRemove("--debug-mem"); + Int_t kDebugEnv = arg.HasOnlyAndRemove("--debug-env") ? 1 : 0; + kDebugEnv = arg.GetIntAndRemove("--debug-env=", kDebugEnv); + + const Bool_t kQuit = arg.HasOnlyAndRemove("-q"); + const Bool_t kBatch = arg.HasOnlyAndRemove("-b"); + const Bool_t kOverwrite = arg.HasOnlyAndRemove("-f"); + + const Bool_t kRefill = arg.HasOnlyAndRemove("--refill"); + + if (arg.GetNumOptions()>0) + { + gLog << warn << "WARNING - Unknown commandline options..." << endl; + arg.Print("options"); + gLog << endl; + return 2; + } + + // + // check for the right usage of the program + // + if (arg.GetNumArguments()<2 || arg.GetNumArguments()>3) + { + Usage(); + return 2; + } + + // + // Setup sequence file and check for its existance + // + const TString kInfile = arg.GetArgumentStr(0); + const TString kDataset = arg.GetArgumentStr(1); + const TString kOutfile = arg.GetArgumentStr(2); + + if (gSystem->AccessPathName(kInfile, kFileExists)) + { + gLog << err << "Sorry, sequences file '" << kInfile << "' doesn't exist." << endl; + return 2; + } + if (gSystem->AccessPathName(kDataset, kFileExists)) + { + gLog << err << "Sorry, dataset file '" << kDataset << "' doesn't exist." << endl; + return 2; + } + + if (gSystem->AccessPathName(kOutfile, kWritePermission)) + { + gLog << err << "Sorry, no write permission for '" << kOutfile << "' doesn't exist." << endl; + return 2; + } + + if (kDebugMem) + TObject::SetObjectStat(kTRUE); + + // + // Setup sequence and check its validity + // + MDataSet seq(kDataset); + if (kPrintSeq || kPrintFiles) + { + gLog << all; + gLog.Separator(kDataset); + seq.Print(kPrintFiles?"files":""); + gLog << endl; + } + if (!seq.IsValid()) + { + gLog << err << "Sequences read but not valid!" << endl << endl; + return 2; + } + + // + // Initialize root + // + MArray::Class()->IgnoreTObjectStreamer(); + MParContainer::Class()->IgnoreTObjectStreamer(); + + TApplication app("Sponde", &argc, argv); + if (!gROOT->IsBatch() && !gClient || gROOT->IsBatch() && !kBatch) + { + gLog << err << "Bombing... maybe your DISPLAY variable is not set correctly!" << endl; + return 1; + } + + // + // Update frequency by default = 1Hz + // + MStatusDisplay *d = new MStatusDisplay; + + // From now on each 'Exit' means: Terminate the application + d->SetBit(MStatusDisplay::kExitLoopOnExit); + d->SetTitle(kDataset); + + // + // Calculate pedestal for pedestal-calculation and calibration + // + MEnv env(kConfig); + + MJSpectrum job(Form("MJSpectrum #%d", seq.GetNumAnalysis())); + job.SetEnv(&env); + job.SetEnvDebug(kDebugEnv); + job.SetDisplay(d);; + job.SetOverwrite(kOverwrite); + job.SetPathOut(kOutfile); + job.SetPathIn(kInfile); + job.EnableRefilling(kRefill); + + if (!job.Process(seq)) + { + gLog << err << "Calculation of spectrum failed." << endl << endl; + return 2; + } + if (kDebugEnv>0) + env.PrintUntouched(); + + if (!job.GetDisplay()) + { + gLog << warn << "Display closed by user... execution aborted." << endl << endl; + return 1; + } + + if (kBatch || kQuit) + delete d; + else + { + // From now on each 'Close' means: Terminate the application + d->SetBit(MStatusDisplay::kExitLoopOnClose); + + // Wait until the user decides to exit the application + app.Run(kFALSE); + } + + if (TObject::GetObjectStat()) + { + TObject::SetObjectStat(kFALSE); + gObjectTable->Print(); + } + + return 0; +}