Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6937) +++ trunk/MagicSoft/Mars/Changelog (revision 6938) @@ -21,4 +21,23 @@ -*-*- END OF LINE -*-*- + + 2005/04/13 Thomas Bretz + + * mbase/MStatusDisplay.[h,cc]: + - implemented new Image types which can be written + + * mhbase/MBinning.[h,cc]: + - added some more constructors which can set the binning directly + + * mhflux/MHEnergyEst.h: + - added getter for fHEnergy + + * mhflux/FluxLinkDef.h, mhflux/Makefile: + - added MHCollectionArea + + * mhflux/MHCollectionArea.[h,cc]: + - added + + 2005/04/13 Daniela Dorner Index: trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc (revision 6937) +++ trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc (revision 6938) @@ -31,5 +31,5 @@ // (see Read and Write) or printed as a postscript file (see SaveAsPS). // -// To write gif files of C-Macros use SaveAsGif() or SaveAsC(). +// To write gif files of C-Macros use SaveAsGif()/SaveAsPNG() or SaveAsC(). // Direct printing to the default printer (via lpr) can be done by // PrintToLpr(). @@ -193,5 +193,8 @@ filemenu->AddSeparator(); filemenu->AddEntry("Save status.&ps", kFileSaveAsPS); + filemenu->AddEntry("Save status.&png", kFileSaveAsPNG); filemenu->AddEntry("Save status.&gif", kFileSaveAsGIF); + filemenu->AddEntry("Save status.&jpg", kFileSaveAsJPG); + filemenu->AddEntry("Save status.&xpm", kFileSaveAsXPM); filemenu->AddEntry("Save status.&C", kFileSaveAsC); filemenu->AddEntry("Save status.&root", kFileSaveAsRoot); @@ -217,5 +220,8 @@ tabmenu->AddSeparator(); tabmenu->AddEntry("Save tab-i.&ps", kTabSaveAsPS); + tabmenu->AddEntry("Save tab-i.&png", kTabSaveAsPNG); tabmenu->AddEntry("Save tab-i.&gif", kTabSaveAsGIF); + tabmenu->AddEntry("Save tab-i.&jpg", kTabSaveAsJPG); + tabmenu->AddEntry("Save tab-i.&xpm", kTabSaveAsXPM); tabmenu->AddEntry("Save tab-i.&C", kTabSaveAsC); tabmenu->AddEntry("Save tab-i.&root", kTabSaveAsRoot); @@ -1198,8 +1204,20 @@ return kTRUE; + case kFileSaveAsPNG: + SaveAsPNG(); + return kTRUE; + case kFileSaveAsGIF: SaveAsGIF(); return kTRUE; + case kFileSaveAsXPM: + SaveAsXPM(); + return kTRUE; + + case kFileSaveAsJPG: + SaveAsJPG(); + return kTRUE; + case kFileSaveAsC: SaveAsC(); @@ -1222,6 +1240,18 @@ return kTRUE; + case kTabSaveAsPNG: + SaveAsPNG(fTab->GetCurrent()); + return kTRUE; + case kTabSaveAsGIF: SaveAsGIF(fTab->GetCurrent()); + return kTRUE; + + case kTabSaveAsXPM: + SaveAsXPM(fTab->GetCurrent()); + return kTRUE; + + case kTabSaveAsJPG: + SaveAsJPG(fTab->GetCurrent()); return kTRUE; @@ -1869,5 +1899,12 @@ } -void MStatusDisplay::AddExtension(TString &name, const TString &ext, Int_t num) const +// -------------------------------------------------------------------------- +// +// If the filename name doesn't end with ext, ext is added to the end. +// If name.IsNull() "status" is assumed and the a number (>0) is added +// as "status-6". +// The extension is returned. +// +const TString &MStatusDisplay::AddExtension(TString &name, const TString &ext, Int_t num) const { if (name.IsNull()) @@ -1882,8 +1919,10 @@ if (name.EndsWith("."+ext)) - return; + return ext; name += "."; name += ext; + + return ext; } @@ -2170,12 +2209,14 @@ } -Bool_t MStatusDisplay::SaveAsGIF(Int_t num, TString name) +Bool_t MStatusDisplay::SaveAsImage(Int_t num, TString name, TImage::EImageFileTypes type) { if (gROOT->IsBatch()) { - *fLog << warn << "Sorry, writing gif-files is not available in batch mode." << endl; + *fLog << warn << "Sorry, writing image-files is not available in batch mode." << endl; + //*fLog << warn << "Sorry, writing gif-files is not available in batch mode." << endl; return 0; } - SetStatusLine1("Writing GIF file..."); + //SetStatusLine1("Writing GIF file..."); + SetStatusLine1("Writing image file... "); SetStatusLine2(""); @@ -2186,8 +2227,63 @@ } - AddExtension(name, "gif", num); + TString ext; + switch (type) + { + case TImage::kXpm: + case TImage::kZCompressedXpm: + ext = AddExtension(name, "xpm", num); + break; +// case TImage::kGZCompressedXpm: +// ext = AddExtension(name, "xpm.gz", num); +// break; + case TImage::kPng: + ext = AddExtension(name, "png", num); + break; + case TImage::kJpeg: + ext = AddExtension(name, "jpg", num); + break; +// case TImage::kXcf: +// ext = AddExtension(name, "xcf", num); +// break; +// case TImage::kPpm: +// ext = AddExtension(name, "ppm", num); +// break; +// case TImage::kPnm: +// ext = AddExtension(name, "pnm", num); +// break; +// case TImage::kBmp: +// ext = AddExtension(name, "bmp", num); +// break; +// case TImage::kIco: +// ext = AddExtension(name, "ico", num); +// break; +// case TImage::kCur: +// ext = AddExtension(name, "cur", num); +// break; + case TImage::kGif: + ext = AddExtension(name, "gif", num); + break; +// case TImage::kTiff: +// ext = AddExtension(name, "tif", num); +// break; +// case TImage::kXbm: +// ext = AddExtension(name, "xbm", num); +// break; +// case TImage::kFits: +// ext = AddExtension(name, "fits", num); +// break; +// case TImage::kTga: +// ext = AddExtension(name, "tga", num); +// break; +// case TImage::kXml: +// ext = AddExtension(name, "xml", num); +// break; + default: + *fLog << warn << "Sorry, unknown or unsupported file type..." << endl; + return 0; + } if (num<0) - *fLog << inf << "Writing gif-Files..." << endl; + *fLog << inf << "Writing " << ext << "-Files..." << endl; TPad *padsav = (TPad*)gPad; @@ -2239,10 +2335,13 @@ c->Draw(); - c->SaveAs(writename); - /* - n->Draw(); - n->SaveAs(writename); - delete n; - */ + if (type==TImage::kGif) + c->SaveAs(writename); // FIXME: Seems not to work well in TImage! (root 3.10/02) + else + { + TImage *img = TImage::Create(); + img->FromPad(c); + img->WriteImage(writename, type); + delete img; + } if (num<0) Index: trunk/MagicSoft/Mars/mbase/MStatusDisplay.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MStatusDisplay.h (revision 6937) +++ trunk/MagicSoft/Mars/mbase/MStatusDisplay.h (revision 6938) @@ -12,4 +12,8 @@ #ifndef ROOT_TTimer #include +#endif + +#ifndef ROOT_TImage +#include #endif @@ -38,11 +42,13 @@ // kFile kFileBrowser, kFileCanvas, kFileOpen, kFileSave, kFileSaveAs, kFileSaveAsPS, - kFileSaveAsRoot, kFileSaveAsGIF, kFileSaveAsC, kFilePrint, - kFilePrinterName, kFileClose, kFileExit, kFileReset, + kFileSaveAsRoot, kFileSaveAsPNG, kFileSaveAsGIF, kFileSaveAsJPG, + kFileSaveAsXPM, kFileSaveAsC, kFilePrint, kFilePrinterName, + kFileClose, kFileExit, kFileReset, // kLoop kLoopNone, kLoopStop, // kTab - kTabSave, kTabSaveAs, kTabSaveAsPS, kTabSaveAsRoot, kTabSaveAsGIF, - kTabSaveAsC, kTabPrint, kTabNext, kTabPrevious, kTabRemove, + kTabSave, kTabSaveAs, kTabSaveAsPS, kTabSaveAsRoot, kTabSaveAsPNG, + kTabSaveAsGIF, kTabSaveAsJPG, kTabSaveAsXPM, kTabSaveAsC, + kTabPrint, kTabNext, kTabPrevious, kTabRemove, // kSize kSize640, kSize768, kSize800, kSize960, kSize1024, kSize1152, @@ -124,5 +130,5 @@ Bool_t Display(const TObjArray &list, const char *tab=0); - void AddExtension(TString &name, const TString &ext, Int_t num) const; + const TString &AddExtension(TString &name, const TString &ext, Int_t num) const; void UpdatePSHeader(const TString &name) const; @@ -132,4 +138,6 @@ TRootEmbeddedCanvas *GetEmbeddedCanvas(TGCompositeFrame *cf) const; + + Bool_t SaveAsImage(Int_t num, TString name, TImage::EImageFileTypes type); public: @@ -187,5 +195,8 @@ Int_t SaveAsPS(TString name="", const TString addon="") { return SaveAsPS(-1, name, addon); } + Bool_t SaveAsPNG(TString name="") { return SaveAsPNG(-1, name); } Bool_t SaveAsGIF(TString name="") { return SaveAsGIF(-1, name); } + Bool_t SaveAsXPM(TString name="") { return SaveAsXPM(-1, name); } + Bool_t SaveAsJPG(TString name="") { return SaveAsJPG(-1, name); } Bool_t SaveAsC(TString name="") { return SaveAsC(-1, name); } Int_t SaveAsRoot(TString name="") { return SaveAsRoot(-1, name); } @@ -193,5 +204,8 @@ Int_t SaveAsPS(Int_t num, TString name="", const TString addon=""); - Bool_t SaveAsGIF(Int_t num, TString name=""); + Bool_t SaveAsPNG(Int_t num, TString name="") { return SaveAsImage(num, name, TImage::kPng); } + Bool_t SaveAsGIF(Int_t num, TString name="") { return SaveAsImage(num, name, TImage::kGif); } + Bool_t SaveAsXPM(Int_t num, TString name="") { return SaveAsImage(num, name, TImage::kXpm); } + Bool_t SaveAsJPG(Int_t num, TString name="") { return SaveAsImage(num, name, TImage::kJpeg); } Bool_t SaveAsC(Int_t num, TString name=""); Int_t SaveAsRoot(Int_t num, TString name=""); Index: trunk/MagicSoft/Mars/mhbase/MBinning.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MBinning.cc (revision 6937) +++ trunk/MagicSoft/Mars/mhbase/MBinning.cc (revision 6938) @@ -91,4 +91,33 @@ // -------------------------------------------------------------------------- // +// Initialize Binning from an axis of a TH1. If no title given, +// a title combined from the axis titles and the TH1 title is +// used. +// +MBinning::MBinning(const TH1 &h, const Char_t axis, const char *name, const char *title) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? Form("%s;%s;%s;%s", h.GetTitle(), + h.GetXaxis()->GetTitle(), + h.GetYaxis()->GetTitle(), + h.GetZaxis()->GetTitle()) : gsDefTitle.Data(); + + SetEdges(h, axis); +} + +// -------------------------------------------------------------------------- +// +// Initialize Binning from TAxis. +// +MBinning::MBinning(const TAxis &axis, const char *name, const char *title) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + SetEdges(axis); +} + +// -------------------------------------------------------------------------- +// // Search in the parameter list for the binning with name "name". If found, // set the edges and title accordingly. Default is name of object. @@ -412,8 +441,8 @@ *fLog << all; *fLog << GetDescriptor() << ": nbins=" << GetNumBins() << " ["; - *fLog << GetEdgeLo() << ", " << GetEdgeHi() << "] "; + *fLog << GetEdgeLo() << ", " << GetEdgeHi() << "] <"; switch (fType) { - case kIsDefault: *fLog << "deafult"; break; + case kIsDefault: *fLog << "default"; break; case kIsLinear: *fLog << "linear"; break; case kIsLogarithmic: *fLog << "logarithmic"; break; Index: trunk/MagicSoft/Mars/mhbase/MBinning.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MBinning.h (revision 6937) +++ trunk/MagicSoft/Mars/mhbase/MBinning.h (revision 6938) @@ -41,4 +41,6 @@ MBinning(Int_t nbins, Axis_t lo, Axis_t hi, const char *name=0, const char *opt="", const char *title=NULL); MBinning(const MBinning &bins) { SetEdges(bins); } + MBinning(const TH1 &h, const Char_t axis='x', const char *name=0, const char *title=0); + MBinning(const TAxis &a, const char *name=0, const char *title=0); void Copy(TObject &named) const Index: trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h (revision 6937) +++ trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h (revision 6938) @@ -11,4 +11,5 @@ #pragma link C++ class MHFalseSource+; #pragma link C++ class MHEffectiveOnTime+; +#pragma link C++ class MHCollectionArea+; #endif Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 6938) +++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 6938) @@ -0,0 +1,396 @@ +/* ======================================================================== *\ +! +! * +! * 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 12/2000 +! Author(s): Harald Kornmayer 1/2001 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHCollectionArea +// +////////////////////////////////////////////////////////////////////////////// +#include "MHCollectionArea.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MH.h" +#include "MBinning.h" + +#include "MMcEvt.hxx" +#include "MMcRunHeader.hxx" +#include "MMcCorsikaRunHeader.h" + +#include "MParList.h" +#include "MParameters.h" + +ClassImp(MHCollectionArea); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Creates the three necessary histograms: +// - selected showers (input) +// - all showers (input) +// - collection area (result) +// +MHCollectionArea::MHCollectionArea(const char *name, const char *title) + : fMcEvt(0), fEnergy(0), fMcAreaRadius(300.), fIsExtern(kFALSE) +{ + // initialize the histogram for the distribution r vs E + // + // we set the energy range from 2 Gev to 20000 GeV (in log 4 orders + // of magnitude) and for each order we take 25 subdivision --> 100 xbins + // + // we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins + // + fName = name ? name : "MHCollectionArea"; + fTitle = title ? title : "Collection Area vs. Energy/Theta"; + + fHistSel.SetName("SelEvts"); + fHistSel.SetTitle("Number of Events after cuts"); + fHistSel.SetXTitle("\\Theta [deg]"); + fHistSel.SetYTitle("E [GeV]"); + fHistSel.SetDirectory(NULL); + fHistSel.UseCurrentStyle(); + + fHistAll.SetName("AllEvts"); + fHistAll.SetTitle("Number of events produced"); + fHistAll.SetXTitle("\\Theta [deg]"); + fHistAll.SetYTitle("E [GeV]"); + fHistAll.SetDirectory(NULL); + fHistAll.UseCurrentStyle(); + + fHEnergy.SetName("CollEnergy"); + fHEnergy.SetTitle("Collection Area vs Energy"); + fHEnergy.SetXTitle("E [GeV]"); + fHEnergy.SetYTitle("A [m^{2}]"); + fHEnergy.SetDirectory(NULL); + fHEnergy.UseCurrentStyle(); + + MBinning binsa, binse, binst; + binse.SetEdgesLog(15, 10, 1000000); + binst.SetEdgesCos(50, 0, 60); + + binse.Apply(fHEnergy); + + MH::SetBinning(&fHistSel, &binst, &binse); + MH::SetBinning(&fHistAll, &binst, &binse); +} + +// -------------------------------------------------------------------------- +// +// Calculate the Efficiency (collection area) and set the 'ReadyToSave' +// flag +// +void MHCollectionArea::CalcEfficiency() +{ + TH1D *hsel = fHistSel.ProjectionY(); + TH1D *hall = fHistAll.ProjectionY(); + + const Int_t nbinx = hsel->GetNbinsX(); + + // ----------------------------------------------------------- + // + // Impact parameter range: TO BE FIXED! Impact parameter shoud be + // read from run header, but it is not yet in!! + // + const Float_t r1 = 0; + const Float_t r2 = fMcAreaRadius; + + //*fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of " << fMCAreaRadius << " meters for the MC events. Check that this is the true one!" <GetBinContent(ix); + + if (Na <= 0) + continue; + + const Float_t Ns = hsel->GetBinContent(ix); + + // Since Na is an estimate of the total number of showers generated + // in the energy bin, it may happen that Ns (triggered showers) is + // larger than Na. In that case, the bin is skipped: + + if (Na < Ns) + continue; + + const Double_t eff = Ns/Na; + + const Double_t efferr = TMath::Sqrt((1.-eff)*Ns)/Na; + + const Float_t col_area = eff * total_area; + const Float_t col_area_error = efferr * total_area; + + fHEnergy.SetBinContent(ix, col_area); + fHEnergy.SetBinError(ix, col_area_error); + } + + delete hsel; + delete hall; +} + +Bool_t MHCollectionArea::SetupFill(const MParList *pl) +{ + fHistSel.Reset(); + if (!fIsExtern) + fHistAll.Reset(); + + fMcEvt = (MMcEvt*)pl->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << "MMcEvt not found... abort." << endl; + return kFALSE; + } + + fEnergy = (MParameterD*)pl->FindObject("MEnergyEst", "MParameterD"); + if (!fEnergy) + { + *fLog << err << "MEnergyEst [MParameterD] not found... abort." << endl; + return kFALSE; + } + + MBinning binst, binse; + binst.SetEdges(fHistAll, 'x'); + binse.SetEdges(fHistAll, 'y'); + + //if (!fIsExtern) + { + MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning"); + if (bins) + binst.SetEdges(*bins); + + bins = (MBinning*)pl->FindObject("BinningEnergy", "MBinning"); + if (bins) + binse.SetEdges(*bins); + } + + binse.Apply(fHEnergy); + + MH::SetBinning(&fHistSel, &binst, &binse); + MH::SetBinning(&fHistAll, &binst, &binse); + + return kTRUE; +} + +Bool_t MHCollectionArea::ReInit(MParList *plist) +{ + if (fIsExtern) + return kTRUE; + + MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader"); + if (!runheader) + { + *fLog << err << "MMcRunHeader not found... abort." << endl; + return kFALSE; + } + + fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers(); + *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl; + + /* + if (fTheta>=0 && fTheta!=runheader->GetTelesTheta()) + { + *fLog << warn << "Warning - Read files have different TelesTheta ("; + *fLog << fTheta << ", " << runheader->GetTelesTheta() << ")..." << endl; + } + fTheta = runheader->GetTelesTheta(); + */ + if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion()) + *fLog << warn << "Warning - Read files have different Corsika versions..." << endl; + fCorsikaVersion = runheader->GetCorsikaVersion(); + + fAllEvtsTriggered |= runheader->GetAllEvtsTriggered(); + *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl; + + MMcCorsikaRunHeader *crh = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader"); + if (!crh) + { + *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl; + return kFALSE; + } + + // + // Calculate approximately the original number of events in each + // energy bin: + // + const Float_t emin = crh->GetELowLim(); + const Float_t emax = crh->GetEUppLim(); + const Float_t expo = 1 + crh->GetSlopeSpec(); + const Float_t k = runheader->GetNumSimulatedShowers() / + (pow(emax,expo) - pow(emin,expo)); + + const Int_t nbiny = fHistAll.GetNbinsY(); + + TAxis &axe = *fHistAll.GetYaxis(); + for (Int_t i = 1; i <= nbiny; i++) + { + const Float_t e1 = axe.GetBinLowEdge(i); + const Float_t e2 = axe.GetBinLowEdge(i+1); + + if (e1 < emin || e2 > emax) + continue; + + const Float_t events = k * (pow(e2, expo) - pow(e1, expo)); + // + // We fill the i-th energy bin, with the total number of events + // Second argument of Fill would be impact parameter of each + // event, but we don't really need it for the collection area, + // so we just put a dummy value (1.) + // + + const Float_t energy = (e1+e2)/2.; + fHistAll.Fill(20, energy, events); + // you have MMcRunHeader.fShowerThetaMin and MMcRunHeader.fShowerThetaMax + } + + return kTRUE; +} + +void MHCollectionArea::Paint(Option_t *option) +{ + TVirtualPad *pad = gPad; + + TH1 *h=0; + + pad->cd(1); + //if (gPad->FindObject("ProjAllX")) + // fHistAll.ProjectionX("ProjAllX", -1, 9999, "E"); + if (gPad->FindObject("ProjSelX")) + fHistSel.ProjectionX("ProjSelX", -1, 9999, "E"); + + pad->cd(2); + if (gPad->FindObject("ProjAllY")) + h=fHistAll.ProjectionY("ProjAllY", -1, 9999, "E"); + if (gPad->FindObject("ProjSelY")) + fHistSel.ProjectionY("ProjSelY", -1, 9999, "E"); + + if (h->GetMaximum()>0) + { + gPad->SetLogx(); + gPad->SetLogy(); + } + pad->cd(4); + CalcEfficiency(); + if (fHEnergy.GetMaximum()>0) + { + gPad->SetLogx(); + gPad->SetLogy(); + } +} + +void MHCollectionArea::Draw(Option_t *option) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + // Do the projection before painting the histograms into + // the individual pads + AppendPad(); + + pad->SetBorderMode(0); + pad->Divide(2,2); + + TH1D *h=0; + + if (fHistSel.GetNbinsX()>1) + { + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + /* + h = fHistAll.ProjectionX("ProjAllX", -1, 9999, "E"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetDirectory(NULL); + h->SetLineColor(kGreen); + h->SetBit(kCanDelete); + h->Draw(); + */ + h = fHistSel.ProjectionX("ProjSelX", -1, 9999, "E"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetDirectory(NULL); + h->SetLineColor(kRed); + h->SetBit(kCanDelete); + h->Draw("hist"/*"same"*/); + } + else + delete pad->GetPad(1); + + if (fHistSel.GetNbinsY()>1) + { + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + + h = fHistAll.ProjectionY("ProjAllY", -1, 9999, "E"); + h->SetDirectory(NULL); + h->SetLineColor(kGreen); + h->SetXTitle("E [GeV]"); + h->SetBit(kCanDelete); + h->Draw(); + + h = fHistSel.ProjectionY("ProjSelY", -1, 9999, "E"); + h->SetDirectory(NULL); + h->SetLineColor(kRed); + h->SetBit(kCanDelete); + h->Draw("same"); + } + else + delete pad->GetPad(2); + + if (fHEnergy.GetNbinsX()>1) + { + pad->cd(4); + gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + fHEnergy.Draw(); + } + else + delete pad->GetPad(4); +/* + if (fHistAll.GetNbinsY()>1) + { + pad->cd(4); + } + else + delete pad->GetPad(4); + */ +} + +Bool_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight) +{ + const Double_t energy = fMcEvt->GetEnergy()/*fEnergy->GetVal()*/; + const Double_t theta = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg(); + + //*fLog << energy << " " << theta << endl; + + fHistSel.Fill(theta, energy); + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h (revision 6938) +++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h (revision 6938) @@ -0,0 +1,76 @@ +#ifndef MARS_MHCollectionArea +#define MARS_MHCollectionArea + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH1 +#include +#endif + +#ifndef ROOT_TH2 +#include +#endif + +class MMcEvt; +class MParameterD; + +class MHCollectionArea : public MH +{ +private: + MMcEvt *fMcEvt; + MParameterD *fEnergy; + + TH2D fHistSel; + TH2D fHistAll; + + TH1D fHEnergy; + + Float_t fMcAreaRadius; // [m] Radius of circle within which the cores of Corsika events have been uniformly distributed. + + UInt_t fTotalNumSimulatedShowers; + UInt_t fCorsikaVersion; + Bool_t fAllEvtsTriggered; + Bool_t fIsExtern; + + void Calc(TH2D &hsel, TH2D &hall); + +public: + MHCollectionArea(const char *name=NULL, const char *title=NULL); + + Bool_t SetupFill(const MParList *pList); + Bool_t ReInit(MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t weight=1); + + void Draw(Option_t *option=""); + void Paint(Option_t *option=""); + + void CalcEfficiency(); + + void SetHistAll(const TH2D &h) { h.Copy(fHistAll); fIsExtern=kTRUE; } + + TH2D &GetHistAll() { return fHistAll; } + const TH2D &GetHistAll() const { return fHistAll; } + +/* + void DrawAll(Option_t *option=""); + void DrawSel(Option_t *option=""); + + const TH1D *GetHist() { return fHistCol; } + const TH1D *GetHist() const { return fHistCol; } + + + + void CalcEfficiency(); + void CalcEfficiency2(); + + void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall); + void Calc(const MHMcEfficiency &heff); + */ + void SetMcAreaRadius(Float_t x) { fMcAreaRadius = x; } + + ClassDef(MHCollectionArea, 1) // Data Container to calculate Collection Area +}; + +#endif Index: trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h (revision 6937) +++ trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h (revision 6938) @@ -48,4 +48,6 @@ MHEnergyEst(const char *name=NULL, const char *title=NULL); + const TH3D &GetHEnergy() const { return fHEnergy; } + void InitMapping(MHMatrix *mat); void StopMapping(); Index: trunk/MagicSoft/Mars/mhflux/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhflux/Makefile (revision 6937) +++ trunk/MagicSoft/Mars/mhflux/Makefile (revision 6938) @@ -29,4 +29,5 @@ MHEnergyEst.cc \ MHEffectiveOnTime.cc \ + MHCollectionArea.cc \ MHFalseSource.cc