Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6489) +++ trunk/MagicSoft/Mars/Changelog (revision 6491) @@ -21,4 +21,19 @@ -*-*- END OF LINE -*-*- + + 2005/02/15 Abelardo Moralejo + + * mmontecarlo/MMcCollectionAreaCalc.[h,cc] + * mhistmc/MHMcCollectionArea.[h,cc] + - Changed the way of calculating final effective area for data + analysis. The new approach requires the use of MC files produced + with the current CVS version of camera. We now make use of the + true total number of produced MC events, and allow for the + setting of a "tentative" differential gamma spectrum to be used + in the calculation of effective areas. + + * macros/collarea.C + - Adapted to the new way of calculating effective areas. + 2005/02/15 Thomas Bretz Index: trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc (revision 6489) +++ trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc (revision 6491) @@ -1,26 +1,27 @@ /* ======================================================================== *\ -! -! * -! * 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 -! -! -\* ======================================================================== */ + ! + ! * + ! * 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 + ! Author(s): Abelardo Moralejo 2/2005 + ! + ! Copyright: MAGIC Software Development, 2000-2005 + ! + ! + \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// @@ -33,76 +34,113 @@ #include +#include #include +#include +#include +#include #include "MH.h" #include "MBinning.h" -#include "MHMcEfficiency.h" -#include "MHMcEnergyImpact.h" #include "MLog.h" #include "MLogManip.h" + ClassImp(MHMcCollectionArea); using namespace std; -// -------------------------------------------------------------------------- -// -// Creates the three necessary histograms: +//////////////////////////////////////////////////////////////////////////////// +// +// Constructor. Creates the three necessary histograms: // - selected showers (input) // - all showers (input) // - collection area (result) // -MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title) - : fMCAreaRadius(300.) +MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title): + fImpactBins(50), fImpactMax(500.), fMinEvents(10) { - // 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 : "MHMcCollectionArea"; - fTitle = title ? title : "Collection Area vs. Energy"; - - MBinning binsx; - MBinning binsy; - binsx.SetEdgesLog(100, 2., 20000); - binsy.SetEdges (50, 0, 500); - - fHistAll = new TH2D; - fHistSel = new TH2D; - fHistCol = new TH1D; + fName = name ? name : "CollectionArea"; + fTitle = title ? title : "Collection Area vs. Theta vs. Energy"; + + // + // Initialize the histogram for the distribution + // Theta vs impact parameter vs E (z, y, x) + // + // As default we set the energy range from 2 Gev to 20000 GeV (in log 4 + // orders of magnitude) and for each order we take 25 subdivisions --> + // 100 xbins. We set the radius range from 0 m to 500 m with 10 m bin --> + // 50 ybins. We make bins equally spaced in cos(theta) + // + // The coarse binning (of fHistColCoarse) is not set by default, the + // PreProcess of mmc/MMcCollectionAreaCalc will do it with the binnings + // found in the parameter list. + // + + MBinning binsx; + MBinning binsy; + MBinning binsz; + + Int_t nbins = 32; + TArrayD edges(nbins+1); + + edges[0] = 0; + + for(int i = 0; i < nbins; i++) + { + Double_t x = 1 - i*0.01; // x = cos(theta) + edges[i+1] = acos(x-0.005)*kRad2Deg; + } + + binsx.SetEdgesLog(100, 2., 20000); // Energy [GeV] + binsy.SetEdges (fImpactBins, 0, fImpactMax); // Impact parameter [m] + binsz.SetEdges (edges); // Theta [deg] + + fHistAll = new TH3D(); + fHistSel = new TH3D(); + fHistCol = new TH2D(); + fHistColCoarse = new TH2D(); - MH::SetBinning(fHistAll, &binsx, &binsy); - MH::SetBinning(fHistSel, &binsx, &binsy); - - fHistCol->SetName(fName); - fHistAll->SetName("AllEvents"); - fHistSel->SetName("SelectedEvents"); - - fHistCol->SetTitle(fTitle); - fHistAll->SetTitle("All showers - Radius vs Energy distribution"); - fHistSel->SetTitle("Selected showers - Radius vs Energy distribution"); - - fHistAll->SetDirectory(NULL); - fHistSel->SetDirectory(NULL); - fHistCol->SetDirectory(NULL); - - fHistAll->UseCurrentStyle(); - fHistSel->UseCurrentStyle(); - fHistCol->UseCurrentStyle(); - - fHistAll->SetXTitle("E [GeV]"); - fHistAll->SetYTitle("r [m]"); - fHistAll->SetZTitle("N"); - - fHistSel->SetXTitle("E [GeV]"); - fHistSel->SetYTitle("r [m]"); - fHistSel->SetYTitle("N"); - - fHistCol->SetXTitle("E [GeV]"); - fHistCol->SetYTitle("A [m^{2}]"); + MH::SetBinning(fHistAll, &binsx, &binsy, &binsz); + MH::SetBinning(fHistSel, &binsx, &binsy, &binsz); + + fHistColCoarse->SetName(fName); + fHistCol->SetName("CollAreaFineBins"); + fHistAll->SetName("AllEvents"); + fHistSel->SetName("SelectedEvents"); + + fHistAll->Sumw2(); + fHistSel->Sumw2(); + + fHistColCoarse->SetTitle(fTitle); + fHistCol->SetTitle(fTitle); + fHistAll->SetTitle("All showers - Theta vs Radius vs Energy"); + fHistSel->SetTitle("Selected showers - Theta vs Radius vs Energy"); + + fHistAll->SetDirectory(NULL); + fHistSel->SetDirectory(NULL); + fHistCol->SetDirectory(NULL); + fHistColCoarse->SetDirectory(NULL); + + fHistAll->UseCurrentStyle(); + fHistSel->UseCurrentStyle(); + fHistCol->UseCurrentStyle(); + fHistColCoarse->UseCurrentStyle(); + + fHistAll->SetXTitle("E [GeV]"); + fHistAll->SetYTitle("r [m]"); + fHistAll->SetZTitle("\\theta [\\circ]"); + + fHistSel->SetXTitle("E [GeV]"); + fHistSel->SetYTitle("r [m]"); + fHistSel->SetZTitle("\\theta [\\circ]"); + + fHistCol->SetXTitle("E [GeV]"); + fHistCol->SetYTitle("\\theta [\\circ]"); + fHistCol->SetZTitle("A [m^{2}]"); + + fHistColCoarse->SetXTitle("E [GeV]"); + fHistColCoarse->SetYTitle("\\theta [\\circ]"); + fHistColCoarse->SetZTitle("A [m^{2}]"); } @@ -113,7 +151,38 @@ MHMcCollectionArea::~MHMcCollectionArea() { - delete fHistAll; - delete fHistSel; - delete fHistCol; + delete fHistAll; + delete fHistSel; + delete fHistCol; +} + +// -------------------------------------------------------------------------- +// +// Set the (fine) binnings of histograms fHistAll, fHistSel used in the +// calculations. We do not need to change impact parameter binning. +// +void MHMcCollectionArea::SetBinnings(const MBinning &binsEnergy, + const MBinning &binsTheta) +{ + MBinning binsImpact; + binsImpact.SetEdges(fImpactBins, 0., fImpactMax); + + MH::SetBinning(fHistAll, &binsEnergy, &binsImpact, &binsTheta); + MH::SetBinning(fHistSel, &binsEnergy, &binsImpact, &binsTheta); + + fHistAll->Sumw2(); + fHistSel->Sumw2(); +} + + +// -------------------------------------------------------------------------- +// +// Set the binnings of the histogram fHistColCoarse, the effective areas +// in the coarse bins used in the analysis. +// +// +void MHMcCollectionArea::SetCoarseBinnings(const MBinning &binsEnergy, + const MBinning &binsTheta) +{ + MH::SetBinning(fHistColCoarse, &binsEnergy, &binsTheta); } @@ -122,7 +191,7 @@ // Fill data into the histogram which contains all showers // -void MHMcCollectionArea::FillAll(Double_t energy, Double_t radius) -{ - fHistAll->Fill(energy, radius); +void MHMcCollectionArea::FillAll(Double_t energy, Double_t radius, Double_t theta) +{ + fHistAll->Fill(energy, radius, theta); } @@ -131,210 +200,193 @@ // Fill data into the histogram which contains the selected showers // -void MHMcCollectionArea::FillSel(Double_t energy, Double_t radius) -{ - fHistSel->Fill(energy, radius); -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram with all showers -// -void MHMcCollectionArea::DrawAll(Option_t* option) -{ - if (!gPad) - MH::MakeDefCanvas(fHistAll); - - fHistAll->Draw(option); - - gPad->SetLogx(); - - gPad->Modified(); - gPad->Update(); -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram with the selected showers only. -// -void MHMcCollectionArea::DrawSel(Option_t* option) -{ - if (!gPad) - MH::MakeDefCanvas(fHistSel); - - fHistSel->Draw(option); - - gPad->SetLogx(); - - gPad->Modified(); - gPad->Update(); -} - +void MHMcCollectionArea::FillSel(Double_t energy, Double_t radius, Double_t theta) +{ + fHistSel->Fill(energy, radius, theta); +} + +// -------------------------------------------------------------------------- +// +// Draw +// void MHMcCollectionArea::Draw(Option_t* option) { - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - fHistCol->Draw(option); - - pad->SetLogx(); - - pad->Modified(); - pad->Update(); -} - -// -------------------------------------------------------------------------- -// -// Calculate the Efficiency (collection area) and set the 'ReadyToSave' -// flag -// -void MHMcCollectionArea::CalcEfficiency() -{ - Calc(*fHistSel, *fHistAll); -} - -// -------------------------------------------------------------------------- -// -// Calculate the Efficiency (collection area) and set the 'ReadyToSave' -// flag -// -void MHMcCollectionArea::Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall) -{ - Calc((TH2D&)*mcsel.GetHist(), (TH2D&)*mcall.GetHist()); -} - -// -------------------------------------------------------------------------- -// -// Calculate the Efficiency (collection area) and set the 'ReadyToSave' -// flag -// -void MHMcCollectionArea::Calc(TH2D &hsel, TH2D &hall) -{ - MH::SetBinning(fHistCol, hsel.GetXaxis()); - - fHistCol->Sumw2(); - - TH1D &psel = *hsel.ProjectionX(); - TH1D &pall = *hall.ProjectionX(); - - TH1D &pally = *hall.ProjectionY(); - - Double_t max = 0.; - for (Int_t i = hall.GetNbinsY(); i > 0; i--) - if (pally.GetBinContent(i) > 0) - { - max = pally.GetBinLowEdge(i+1); - break; - } - - fHistCol->Divide(&psel, &pall, TMath::Pi()*max*max, 1); - fHistCol->SetEntries(hsel.GetEntries()); - - delete &psel; - delete &pall; - delete &pally; - - SetReadyToSave(); -} - -// -------------------------------------------------------------------------- -// -// Calculate the Efficiency (collection area) and set the 'ReadyToSave' -// flag -// -void MHMcCollectionArea::CalcEfficiency2() -{ - TH1D &histsel = *fHistSel->ProjectionX(); - TH1D &histall = *fHistAll->ProjectionX(); - - TAxis &xaxis = *histsel.GetXaxis(); - MH::SetBinning(fHistCol, &xaxis); - const Int_t nbinx = xaxis.GetNbins(); - - // ----------------------------------------------------------- - // - // 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!" <SetLogx(); + c1->SetLogz(); + c1->SetGridx(); + c1->SetGridy(); + + fHistCol->Draw("lego2"); + + // + // Averagye Aeff + // + TCanvas *c2 = new TCanvas(); + c2->SetLogx(); + c2->SetLogy(); + c2->SetGridx(); + c2->SetGridy(); + + TH1D* harea = fHistCol->ProjectionX(); + harea->Draw("e1"); + + // + // Plot the Aeff for the different theta + // + TCanvas *c3 = new TCanvas(); + c3->SetLogx(); + c3->SetLogy(); + c3->SetGridx(); + c3->SetGridy(); + + TLegend * leg = new TLegend(0.73,0.65,0.89,0.89); + + TAxis* yaxis = fHistCol->GetYaxis(); + const Int_t nbiny = fHistCol->GetYaxis()->GetNbins(); + + THStack* hs = new THStack("aa","aa"); + + hs->Add(harea,"e1"); + leg->AddEntry(harea,"All","l"); + + for(Int_t iy=1; iy<=nbiny; iy++) { - const Float_t Na = histall.GetBinContent(ix); - - if (Na <= 0) - continue; - - const Float_t Ns = histsel.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 = sqrt((1.-eff)*Ns)/Na; - - const Float_t col_area = eff * total_area; - const Float_t col_area_error = efferr * total_area; - - fHistCol->SetBinContent(ix, col_area); - fHistCol->SetBinError(ix, col_area_error); + + TH1D* h1= fHistCol->ProjectionX(Form("%d",iy),iy,iy); + + if(h1->GetEntries()==0) + continue; + + cout <GetEntries() << endl; + + leg->AddEntry(h1,Form("\\theta = %.0f",yaxis->GetBinCenter(iy)),"l"); + h1->SetLineColor(iy); + hs->Add(h1,"e1"); } - - delete &histsel; - - SetReadyToSave(); -} - -// -------------------------------------------------------------------------- -// -// Calculate the Efficiency (collection area) and set the 'ReadyToSave' -// flag -// -void MHMcCollectionArea::Calc(const MHMcEfficiency &heff) -{ - // - // now calculate the Collection area for different - // energies - // - TH2D &h = (TH2D&)*heff.GetHist(); - - MH::SetBinning(fHistCol, h.GetXaxis()); - - const Int_t nbinx = h.GetXaxis()->GetNbins(); - const Int_t nbiny = h.GetYaxis()->GetNbins(); - - for (Int_t ix=1; ix<=nbinx; ix++) - { - Double_t errA = 0; - Double_t colA = 0; - - for (Int_t iy=1; iy<=nbiny; iy++) - { - TAxis *yaxis = h.GetYaxis(); - - const Double_t r1 = yaxis->GetBinLowEdge(iy); - const Double_t r2 = yaxis->GetBinLowEdge(iy+1); - - const Double_t A = TMath::Pi() * (r2*r2 - r1*r1); - - const Double_t eff = h.GetCellContent(ix, iy); - const Double_t efferr = h.GetCellError(ix, iy); - - colA += eff*A; - errA += A*A * efferr*efferr; - } - - fHistCol->SetBinContent(ix, colA); - fHistCol->SetBinError(ix, sqrt(errA)); - } - - SetReadyToSave(); -} + hs->SetMinimum(1); + + hs->Draw("nostack"); + leg->Draw(); + +} + +// -------------------------------------------------------------------------- +// +// Calculate the collection area and set the 'ReadyToSave' flag +// We first calculate the area in fine energy bins, and then do a +// weighted mean to obtain the area in coarse bins. The weights in +// the coarse bins are intended to account for the effect of the +// energy spectrum in the effective area itself. The weights +// are taken from the tentative differential spectrum dN_gam/dE given +// through the function "spectrum". If no such function is supplied, +// then no weights are applied (and hence the spectrum will be as a +// flat spectrum in dN_gam/dE). Of course we have a "generated" MC +// spectrum, but within each fine bin the differences in spectrum +// should not change the result (if bins are fine enough). With no +// supplied tentative spectrum, each fine bin is weighted equally in +// calculating the area in the coarse bin, and so it is like having a +// flat spectrum. +// +// You can run this Calc procedure on an already existing +// MHMcCollectionArea object, as long as it is filled. +// +void MHMcCollectionArea::Calc(TF1 *spectrum) +{ + // Search last impact parameter bin containing events + // FIXME: this should be done independently for each theta angle. + // + TH1D &himpact = *(TH1D*)fHistAll->Project3D("y"); + + Int_t impbin; + for (impbin = himpact.GetNbinsX(); impbin > 0; impbin--) + if (himpact.GetBinContent(impbin)>0) + break; + + Float_t max_radius = himpact.GetBinLowEdge(impbin); + + Float_t total_area = TMath::Pi()*max_radius*max_radius; + + for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++) + for (Int_t iz = 1; iz <= fHistAll->GetNbinsZ(); iz++) + { + fHistAll->SetBinContent(ix, impbin, iz, 0.); + fHistSel->SetBinContent(ix, impbin, iz, 0.); + fHistAll->SetBinError(ix, impbin, iz, 0.); + fHistSel->SetBinError(ix, impbin, iz, 0.); + } + + TH2D &histsel = *(TH2D*)fHistSel->Project3D("zx,e"); + TH2D &histall = *(TH2D*)fHistAll->Project3D("zx,e"); + // "e" option means that errors are computed! + + + TAxis &xaxis = *histsel.GetXaxis(); + TAxis &yaxis = *histsel.GetYaxis(); + MH::SetBinning(fHistCol, &xaxis, &yaxis); + + cout << "Total considered MC area = pi * " << max_radius + << "^2 square meters" << endl; + + fHistCol->Sumw2(); + fHistCol->Divide(&histsel, &histall, total_area, 1., "b"); + + // + // Now get the effective area in the selected coarse bins. Weight + // the values in the small bins according the supplied tentative + // spectrum, if it has been supplied as argument of Calc. + // + + for (Int_t ibin = 1; ibin <= fHistColCoarse->GetNbinsX(); ibin++) + for (Int_t jbin = 1; jbin <= fHistColCoarse->GetNbinsY(); jbin++) + { + Float_t maxenergy = fHistColCoarse->GetXaxis()->GetBinUpEdge(ibin); + Float_t minenergy = fHistColCoarse->GetXaxis()->GetBinLowEdge(ibin); + + Float_t maxtheta = fHistColCoarse->GetYaxis()->GetBinUpEdge(jbin); + Float_t mintheta = fHistColCoarse->GetYaxis()->GetBinLowEdge(jbin); + + // Fine bins ranges covered by the coarse bin ibin, jbin: + Int_t ibin2max = fHistCol->GetXaxis()->FindBin(maxenergy); + Int_t ibin2min = fHistCol->GetXaxis()->FindBin(minenergy); + + Int_t jbin2max = fHistCol->GetYaxis()->FindBin(maxtheta); + Int_t jbin2min = fHistCol->GetYaxis()->FindBin(mintheta); + + Float_t area = 0.; + Float_t errarea = 0.; + Float_t norm = 0; + + for (Int_t ibin2 = ibin2min; ibin2 <= ibin2max; ibin2++) + { + Float_t weight = spectrum? spectrum-> + Eval(fHistCol->GetXaxis()->GetBinCenter(ibin2)) : 1.; + + for (Int_t jbin2 = jbin2min; jbin2 <= jbin2max; jbin2++) + { + // Skip bins with too few produced MC events + if (histall.GetBinContent(ibin2,jbin2) < fMinEvents) + continue; + + area += weight * fHistCol->GetBinContent(ibin2,jbin2); + norm += weight; + errarea += pow(weight * fHistCol-> + GetBinError(ibin2,jbin2), 2.); + } + } + if (norm > 0.) + { + area /= norm; + errarea = sqrt(errarea)/norm; + } + + fHistColCoarse->SetBinContent(ibin, jbin, area); + fHistColCoarse->SetBinError(ibin, jbin, errarea); + } + + SetReadyToSave(); +} + + Index: trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h (revision 6489) +++ trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h (revision 6491) @@ -6,51 +6,56 @@ #endif +#include "MBinning.h" +#include "TF1.h" + class TH1D; class TH2D; +class TH3D; -class MHMcEfficiency; -class MHMcEnergyImpact; +// +// Version 2: +// +// Now all histograms are saved with the object. Added members +// fImpactMax, fImpactBins, fMinEvents +// class MHMcCollectionArea : public MH { private: - TH2D *fHistAll; //! all simulated showers - TH2D *fHistSel; //! the selected showers + TH3D *fHistAll; // all simulated showers + TH3D *fHistSel; // the selected showers + TH2D *fHistCol; // the collection area in fine bins + TH2D *fHistColCoarse; // the coll. area in coarse bins - TH1D *fHistCol; // the collection area + Int_t fImpactBins; // number of bins in i.p. for histograms + Float_t fImpactMax; // [m] Maximum impact parameter for histograms - Float_t fMCAreaRadius; // [m] Radius of circle within which the cores - // of Corsika events have been uniformly - // distributed. + Int_t fMinEvents; + // minimum number of events in each of the "fine bins" of fHistAll + // so that the bin is used in the averaging for the coarse bins. - void Calc(TH2D &hsel, TH2D &hall); + void Calc(TH3D &hsel, TH3D &hall); public: - MHMcCollectionArea(const char *name=NULL, const char *title=NULL); - ~MHMcCollectionArea(); + MHMcCollectionArea(const char *name=NULL, const char *title=NULL); + ~MHMcCollectionArea(); - void FillAll(Double_t energy, Double_t radius); - void FillSel(Double_t energy, Double_t radius); + void FillAll(Double_t energy, Double_t radius, Double_t theta); + void FillSel(Double_t energy, Double_t radius, Double_t theta); - void DrawAll(Option_t *option=""); - void DrawSel(Option_t *option=""); + const TH2D *GetHist() const { return fHistCol; } + const TH2D *GetHistCoarse() const { return fHistColCoarse; } + const TH3D *GetHistAll() const { return fHistAll; } + const TH3D *GetHistSel() const { return fHistSel; } - const TH1D *GetHist() { return fHistCol; } - const TH1D *GetHist() const { return fHistCol; } + void SetBinnings(const MBinning &binsEnergy, const MBinning &binsTheta); + void SetCoarseBinnings(const MBinning &binsEnergy, const MBinning &binsTheta); - TH2D *GetHistAll() { return fHistAll; } + void Draw(Option_t *option=""); - void Draw(Option_t *option=""); + void Calc(TF1 *spectrum); - 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(MHMcCollectionArea, 1) // Data Container to calculate Collection Area -}; + ClassDef(MHMcCollectionArea, 2) // Data Container to store Collection Area + }; #endif Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc (revision 6489) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc (revision 6491) @@ -17,5 +17,6 @@ ! ! Author(s): Thomas Bretz 12/2000 -! Author(s): Harald Kornmayer 1/2001 +! Author(s): Harald Kornmayer 1/2001 +! Author(s): Abelardo Moralejo 2/2005 ! ! Copyright: MAGIC Software Development, 2000-2002 @@ -26,9 +27,5 @@ ////////////////////////////////////////////////////////////////////////////// // -// MHMcCollectionAreaCalc -// -// Remark: The initialization is mainly done in the ReInit function. -// Please make sure, that you don't use MReadTree when processing -// a file. Use a 'ReInit'-calling task like MReadMarsFile +// MMcCollectionAreaCalc // ////////////////////////////////////////////////////////////////////////////// @@ -42,5 +39,6 @@ #include "MMcEvt.hxx" -#include "MMcTrig.hxx" +#include "MMcEvtBasic.h" + #include "MMcRunHeader.hxx" #include "MMcCorsikaRunHeader.h" @@ -52,181 +50,122 @@ using namespace std; -MMcCollectionAreaCalc::MMcCollectionAreaCalc(const char *input, - const char *name, const char *title) +//////////////////////////////////////////////////////////////////////////////// +// +// Constructor +// +MMcCollectionAreaCalc::MMcCollectionAreaCalc(const char *input, const char *name, + const char *title): + fBinsTheta(0), fBinsEnergy(0), fSpectrum(0) { fName = name ? name : "MMcCollectionAreaCalc"; fTitle = title ? title : "Task to calculate the collection area"; - - fObjName = input ? input : "MMcTrig"; - AddToBranchList(Form("%s.fNumFirstLevel", input)); - - AddToBranchList("MMcEvt.fEnergy"); - AddToBranchList("MMcEvt.fImpact"); - - fCorsikaVersion = 0; - fAllEvtsTriggered = kFALSE; - fTotalNumSimulatedShowers = 0; - fThetaMin = -1.; - fThetaMax = -1.; - } +//////////////////////////////////////////////////////////////////////////////// +// +// PreProcess. +// Create MHMcCollectionArea object if necessary. +// Search in parameter list for MBinning objects with binning in theta and E. +// These contain the coarse binning to be used in the analysis. Then search +// for other necessary input containers: +// if MMcEvt is found, it means we are in the loop over the Events tree, +// and so we must fill the histogram MHMcCollectionArea::fHistSel (using +// MHMcCollectionArea::FillSel). If MMcEvt is not found, it means that we are in +// the loop over the "OriginalMC" tree, containing all the original Corsika events +// produced, and hence we must fill the histogram fHistAll through +// MHMcCollectionArea::FillAll. +// Int_t MMcCollectionAreaCalc::PreProcess (MParList *pList) { - // connect the raw data with this task + fCollArea = (MHMcCollectionArea*)pList->FindCreateObj("MHMcCollectionArea"); - fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); - if (!fMcEvt) + // Look for the binnings of the histograms if they have not been already + // found in a previous loop. + + if (!fBinsTheta) { - *fLog << err << "MMcEvt not found... abort." << endl; - return kFALSE; + fBinsTheta = (MBinning*) pList->FindObject("binsTheta"); + if (!fBinsTheta) + { + *fLog << err << "Coarse Theta binning not found... Aborting." + << endl; + return kFALSE; + } + } + + if (!fBinsEnergy) + { + fBinsEnergy = (MBinning*) pList->FindObject("binsEnergy"); + if (!fBinsEnergy) + { + *fLog << err << "Coarse Energy binning not found... Aborting." + << endl; + return kFALSE; + } + } + + fCollArea->SetCoarseBinnings(*fBinsEnergy, *fBinsTheta); + + + // Look for the input containers + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (fMcEvt) + { + *fLog << inf << "MMcEvt found... I will fill MHMcCollectionArea.fHistSel..." << endl; + return kTRUE; + } + + *fLog << inf << "MMcEvt not found... looking for MMcEvtBasic..." << endl; + + + fMcEvtBasic = (MMcEvtBasic*) pList->FindObject("MMcEvtBasic"); + + if (fMcEvtBasic) + { + *fLog << inf << "MMcEvtBasic found... I will fill MHMcCollectionArea.fHistAll..." << endl; + return kTRUE; + } + + *fLog << err << "MMcEvtBasic not found. Aborting..." << endl; + + return kFALSE; +} + +//////////////////////////////////////////////////////////////////////////////// +// +// Process. Depending on whether fMcEvt or fMcEvtBasic are available, fill +// MHMcCollectionArea::fHistAll, else, if fMcEvt is available, fill +// MHMcCollectionArea::fHistSel +// +Int_t MMcCollectionAreaCalc::Process() +{ + Double_t energy = fMcEvt? fMcEvt->GetEnergy() : fMcEvtBasic->GetEnergy(); + Double_t impact = fMcEvt? fMcEvt->GetImpact()/100. : fMcEvtBasic->GetImpact()/100.; // in m + Double_t theta = fMcEvt? fMcEvt->GetTelescopeTheta()*TMath::RadToDeg() : + fMcEvtBasic->GetTelescopeTheta()*TMath::RadToDeg(); // in deg + + if (fMcEvt) + fCollArea->FillSel(energy, impact, theta); + else + fCollArea->FillAll(energy, impact, theta); + + return kTRUE; +} + +/////////////////////////////////////////////////////////////////////////////// +// +// Postprocess. If both fHistAll and fHistSel are already filled, calculate +// effective areas. Else it means we still have to run one more loop. +// +Int_t MMcCollectionAreaCalc::PostProcess() +{ + if ( ((TH2D*)fCollArea->GetHistAll())->GetEntries() > 0 && + ((TH2D*)fCollArea->GetHistSel())->GetEntries() > 0) + { + *fLog << inf << "Calculation Collection Area..." << endl; + fCollArea->Calc(fSpectrum); } return kTRUE; } - -Bool_t MMcCollectionAreaCalc::ReInit(MParList *plist) -{ - fCollArea=0; - - MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader"); - if (!runheader) - { - *fLog << err << "Error - MMcRunHeader not found... abort." << endl; - return kFALSE; - } - - fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers(); - *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl; - - - if ( (fThetaMin >= 0 && fThetaMin != runheader->GetShowerThetaMin()) || - (fThetaMax >= 0 && fThetaMax != runheader->GetShowerThetaMax()) ) - { - *fLog << warn << "Warning - Read files have different Theta ranges ("; - *fLog << "(" << fThetaMin << ", " << fThetaMax << ") vs (" << - runheader->GetShowerThetaMin() << ", " << - runheader->GetShowerThetaMax() << ")..." << endl; - } - fThetaMin = runheader->GetShowerThetaMin(); - fThetaMax = runheader->GetShowerThetaMax(); - - - 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; - - - fCollArea = (MHMcCollectionArea*)plist->FindCreateObj("MHMcCollectionArea"); - if (!fCollArea) - return kFALSE; - - if (!fAllEvtsTriggered) - { - fMcTrig = (MMcTrig*)plist->FindObject(fObjName); - if (!fMcTrig) - { - *fLog << err << fObjName << " [MMcTrig] not found... anort." << endl; - return kFALSE; - } - return kTRUE; - } - - MMcCorsikaRunHeader *corrunheader = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader"); - if (!corrunheader) - { - *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl; - return kFALSE; - } - - // - // Calculate approximately the original number of events in each - // energy bin: - // - - const Float_t emin = corrunheader->GetELowLim(); - const Float_t emax = corrunheader->GetEUppLim(); - const Float_t expo = 1 + corrunheader->GetSlopeSpec(); - const Float_t k = runheader->GetNumSimulatedShowers() / - (pow(emax,expo) - pow(emin,expo)) ; - - TH2 &hall = *fCollArea->GetHistAll(); - - const Int_t nbinx = hall.GetNbinsX(); - - TAxis &axe = *hall.GetXaxis(); - for (Int_t i = 1; i <= nbinx; 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.; - hall.Fill(energy, 1., events); - } - - return kTRUE; -} - -Int_t MMcCollectionAreaCalc::Process() -{ - const Double_t energy = fMcEvt->GetEnergy(); - const Double_t impact = fMcEvt->GetImpact()/100.; - - // - // This happens for camera files created with Camera 0.5 - // - if (TMath::IsNaN(impact)) - { - *fLog << err << dbginf << "ERROR - Impact=NaN (Not a number)... abort." << endl; - return kERROR; - } - - if (!fAllEvtsTriggered) - { - fCollArea->FillAll(energy, impact); - - if (fMcTrig->GetFirstLevel() <= 0) - return kTRUE; - } - - fCollArea->FillSel(energy, impact); - - return kTRUE; -} - -Int_t MMcCollectionAreaCalc::PostProcess() -{ - if (!fCollArea) - return kTRUE; - - // - // do the calculation of the effectiv area - // - *fLog << inf << "Calculation Collection Area..." << endl; - - if (!fAllEvtsTriggered) - fCollArea->CalcEfficiency(); - else - { - *fLog << inf << "Total number of showers: " << fTotalNumSimulatedShowers << endl; - fCollArea->CalcEfficiency2(); - } - - return kTRUE; -} - Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h (revision 6489) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h (revision 6491) @@ -7,15 +7,28 @@ #include +#include class MParList; class MMcEvt; +class MMcEvtBasic; class MMcTrig; class MHMcCollectionArea; +class MBinning; class MMcCollectionAreaCalc : public MTask { private: - const MMcEvt *fMcEvt; - const MMcTrig *fMcTrig; + const MMcEvt *fMcEvt; + const MMcEvtBasic *fMcEvtBasic; + const MMcTrig *fMcTrig; + + MBinning *fBinsTheta; + MBinning *fBinsEnergy; + // Coarse zenith angle and energy bins used in the analysis + + TF1 *fSpectrum; + // Tentative energy spectrum. This modifies slightly the calculation + // of the effective area (see MHMcCollectionArea::Calc) + MHMcCollectionArea *fCollArea; @@ -23,19 +36,13 @@ TString fObjName; - UInt_t fTotalNumSimulatedShowers; - UInt_t fCorsikaVersion; - Bool_t fAllEvtsTriggered; - Float_t fThetaMin; - Float_t fThetaMax; - - Bool_t ReInit(MParList *plist); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); public: - MMcCollectionAreaCalc(const char *input=NULL, - const char *name=NULL, const char *title=NULL); + MMcCollectionAreaCalc(const char *input = NULL, + const char *name = NULL, const char *title = NULL); + + void SetSpectrum(TF1 *f) { fSpectrum = f; } ClassDef(MMcCollectionAreaCalc, 0) // Task to calculate the collection area histogram