Index: trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc (revision 853) +++ trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc (revision 853) @@ -0,0 +1,150 @@ +/* ======================================================================== *\ +! +! * +! * 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): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +#include "MHMcCollectionArea.h" + +#include +#include + +ClassImp(MHMcCollectionArea) + +MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title) +{ + // + // default constructor + // + + // initialize the histogram for the distribution r vs E + // + // we set the energy range from 1 Gev to 10000 GeV (in log 5 orders + // of magnitude) and for each order we take 10 subdivision --> 50 xbins + // + // we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins + + + *fName = name ? name : "MHMcCollectionArea"; + *fTitle = title ? title : "Data to Calculate Coll-Area"; + + + fHistAll = new TH2D("collAreaAll", "all showers - Radius vs log(E) distribution", + 50, 0., 5., + 50, 0., 500.); + + fHistSel = new TH2D("collAreaSel", "selected showers - Radius vs log(E) distribution", + 50, 0., 5., + 50, 0., 500.); + + fHistCol = new TH1D("collArea", "Collection Area", + 50, 0., 5.); +} + +MHMcCollectionArea::~MHMcCollectionArea() +{ + delete fHistAll; + delete fHistSel; + delete fHistCol; +} + +void MHMcCollectionArea::FillAll(Float_t log10E, Float_t radius) +{ + fHistAll->Fill(log10E, radius); +} + +void MHMcCollectionArea::FillSel(Float_t log10E, Float_t radius) +{ + fHistSel->Fill(log10E, radius); +} + +void MHMcCollectionArea::DrawAll() +{ + fHistAll->Draw(); +} + +void MHMcCollectionArea::DrawSel() +{ + fHistSel->Draw(); +} + +void MHMcCollectionArea::Draw(Option_t* option) +{ + fHistCol->Draw(option); +} + +void MHMcCollectionArea::CalcEfficiency() +{ + // Description! + + // + // first of all calculate the efficency + // do it here by hand to get the right error of efficency + // + const Int_t nbinx = fHistSel->GetXaxis()->GetNbins(); + const Int_t nbiny = fHistSel->GetYaxis()->GetNbins(); + + for (Int_t ix=1; ix<=nbinx; ix++) + { + for (Int_t iy=1; iy<=nbiny; iy++) + { + const Float_t Ns = fHistSel->GetCellContent(ix, iy); + const Float_t Na = fHistAll->GetCellContent(ix, iy); + + if (Na <= 0) + continue; + + const Double_t eff = Ns/Na; + const Double_t err = sqrt(Na + Na*Ns - Ns*Ns - Ns) / (Na*Na); + + fHistSel->SetCellContent(ix, iy, eff); + fHistSel->SetCellError(ix, iy, err); + } + } + + // + // now calculate the Collection area for different + // energies + // + for (Int_t ix=1; ix<=nbinx; ix++) + { + Double_t errA = 0; + Double_t colA = 0; + + for (Int_t iy=1; iy<=nbiny; iy++) + { + const Double_t r1 = fHistSel->GetYaxis()->GetBinLowEdge(iy); + const Double_t r2 = fHistSel->GetYaxis()->GetBinLowEdge(iy+1); + + const Double_t A = kPI * (r2*r2 - r1*r1); + + const Double_t eff = fHistSel->GetCellContent(ix, iy); + const Double_t err = fHistSel->GetCellError(ix, iy); + + colA += eff*A; + errA += A*A * err*err; + } + + fHistCol->SetBinContent(ix, colA); + fHistCol->SetBinError(ix, sqrt(errA)); + } +} Index: trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h (revision 853) +++ trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h (revision 853) @@ -0,0 +1,42 @@ +#ifndef MHMCCOLLECTIONAREA_H +#define MHMCCOLLECTIONAREA_H + +#ifndef MAGIC_H +#include "MAGIC.h" +#endif +#ifndef MPARCONTAINER_H +#include "MParContainer.h" +#endif + +// +// because of some strange reason this cannot be put into MonteCarloIncl +// +#ifndef ROOT_TH1 +#include +#endif + +class TH2D; + +class MHMcCollectionArea : public MParContainer +{ +private: + TH2D *fHistAll ; //! all simulated showers + TH2D *fHistSel ; //! the selected showers + TH1D *fHistCol ; // the collection area + +public: + + MHMcCollectionArea(const char *name=NULL, const char *title=NULL) ; + ~MHMcCollectionArea() ; + + void FillAll(Float_t log10E, Float_t radius) ; + void FillSel(Float_t log10E, Float_t radius) ; + void DrawAll() ; + void DrawSel() ; + void Draw(Option_t* option = "") ; + void CalcEfficiency() ; + + ClassDef(MHMcCollectionArea, 1) // Data Container to calculate Collection Area +}; + +#endif Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc (revision 853) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc (revision 853) @@ -0,0 +1,94 @@ +/* ======================================================================== *\ +! +! * +! * 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): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +#include "MMcCollectionAreaCalc.h" + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MMcEvt.hxx" +#include "MMcTrig.hxx" + +#include "MHMcCollectionArea.h" + +ClassImp(MMcCollectionAreaCalc) + +MMcCollectionAreaCalc::MMcCollectionAreaCalc (const char *name, const char *title) +{ + *fName = name ? name : "MMcCollectionAreaCalc"; + *fTitle = title ? title : "Task to calc the collection area "; +} + +Bool_t MMcCollectionAreaCalc::PreProcess (MParList *pList) +{ + // connect the raw data with this task + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << dbginf << "MMcEvt not found... exit." << endl; + return kFALSE; + } + + fMcTrig = (MMcTrig*)pList->FindObject("MMcTrig"); + if (!fMcTrig) + { + *fLog << dbginf << "MMcTrig not found... exit." << endl; + return kFALSE; + } + + fCollArea = (MHMcCollectionArea*)pList->FindCreateObj("MHMcCollectionArea"); + if (!fCollArea) + return kFALSE; + + return kTRUE; +} + +Bool_t MMcCollectionAreaCalc::Process () +{ + const Float_t energy = log10(fMcEvt->GetEnergy()); + const Float_t impact = fMcEvt->GetImpact()/100.; + + fCollArea->FillAll(energy, impact); + + if (fMcTrig->GetFirstLevel() <= 0) + return kTRUE; + + fCollArea->FillSel(energy, impact); + + return kTRUE; +} + +Bool_t MMcCollectionAreaCalc::PostProcess () +{ + // + // do the calculation of the effectiv area + // + fCollArea->CalcEfficiency(); + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h (revision 853) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h (revision 853) @@ -0,0 +1,32 @@ +#ifndef MMCCOLLECTIONAREACALC_H +#define MMCCOLLECTIONAREACALC_H + +#ifndef MTASK_H +#include "MTask.h" +#endif + +class MParList; +class MMcEvt; +class MMcTrig; +class MHMcCollectionArea; + +class MMcCollectionAreaCalc : public MTask +{ +private: + MMcEvt *fMcEvt; + MMcTrig *fMcTrig; + + MHMcCollectionArea *fCollArea; + +public: + MMcCollectionAreaCalc(const char *name=NULL, const char *title=NULL); + + Bool_t PreProcess(MParList *pList); + Bool_t Process() ; + Bool_t PostProcess() ; + + ClassDef(MMcCollectionAreaCalc, 0) // Task to calculate the collection area histogram +}; + +#endif +