Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 1820) +++ trunk/MagicSoft/Mars/Changelog (revision 1821) @@ -1,3 +1,14 @@ -*-*- END OF LINE -*-*- + + 2003/03/13: Abelardo moralejo + + * mhist/MHMcCT1CollectionArea.[h,cc] + - Added for calculations of collection area for CT1. + + * mmontecarlo/MMcCT1CollectionAreaCalc.[h,cc] + - Added for the same reason. + + * macros/CT1collarea.C + - Uses the above classes 2003/03/12: Abelardo Moralejo Index: trunk/MagicSoft/Mars/mhist/HistLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 1820) +++ trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 1821) @@ -52,4 +52,5 @@ #pragma link C++ class MHMcCollectionArea+; #pragma link C++ class MHMcEnergyMigration+; +#pragma link C++ class MHMcCT1CollectionArea+; #endif Index: trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc (revision 1821) +++ trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc (revision 1821) @@ -0,0 +1,346 @@ +/* ======================================================================== *\ +! +! * +! * 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): A. Moralejo 3/2003 moralejo@pd.infn.it +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHMcCT1CollectionArea // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHMcCT1CollectionArea.h" + +#include +#include + +#include "MH.h" +#include "MBinning.h" + +ClassImp(MHMcCT1CollectionArea); + +// -------------------------------------------------------------------------- +// +// Creates the three necessary histograms: +// - selected showers (input) +// - all showers (input) +// - collection area (result) +// +MHMcCT1CollectionArea::MHMcCT1CollectionArea(const char *name, const char *title) +{ + // initialize the histogram for the distribution r vs E + // + // we set the energy range from 10 Gev to 30000 GeV (in log 4.5 orders + // of magnitude) and for each order we take 10 subdivision --> 45 xbins + // + // we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins + // + fName = name ? name : "MHMcCT1CollectionArea"; + fTitle = title ? title : "Collection Area vs. Energy"; + + MBinning binsx; + MBinning binsy; + binsx.SetEdgesLog(45, 10, 30000); + binsy.SetEdges (50, 0, 500); + + fHistAll = new TH2D; + fHistSel = new TH2D; + fHistCol = new TH1D; + + 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->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}]"); +} + +// -------------------------------------------------------------------------- +// +// Delete the three histograms +// +MHMcCT1CollectionArea::~MHMcCT1CollectionArea() +{ + delete fHistAll; + delete fHistSel; + delete fHistCol; +} + +// -------------------------------------------------------------------------- +// +// Fill data into the histogram which contains the selected showers +// +void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t radius) +{ + fHistSel->Fill(energy, radius); +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram with all showers +// +void MHMcCT1CollectionArea::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 MHMcCT1CollectionArea::DrawSel(Option_t* option) +{ + if (!gPad) + MH::MakeDefCanvas(fHistSel); + + fHistSel->Draw(option); + + gPad->SetLogx(); + + gPad->Modified(); + gPad->Update(); +} + +// -------------------------------------------------------------------------- +// +// Creates a new canvas and draws the histogram into it. +// Be careful: The histogram belongs to this object and won't get deleted +// together with the canvas. +// +TObject *MHMcCT1CollectionArea::DrawClone(Option_t* option) const +{ + TCanvas *c = MH::MakeDefCanvas(fHistCol); + + // + // This is necessary to get the expected bahviour of DrawClone + // + gROOT->SetSelectedPad(NULL); + + fHistCol->DrawCopy(option); + + gPad->SetLogx(); + + c->Modified(); + c->Update(); + + return c; +} + +void MHMcCT1CollectionArea::Draw(Option_t* option) +{ + if (!gPad) + MH::MakeDefCanvas(fHistCol); + + fHistCol->Draw(option); + + gPad->SetLogx(); + + gPad->Modified(); + gPad->Update(); +} + +// +// Calculate the Efficiency (collection area) for the CT1 MC sample +// and set the 'ReadyToSave' flag +// +void MHMcCT1CollectionArea::CalcEfficiency(Float_t theta) +{ + // + // Here we estimate the total number of showers in each energy bin + // from the known the energy range and spectral index of the generated + // showers. This procedure is intended for the CT1 MC files. The total + // number of generated events, collection area, spectral index etc will be + // set here by hand, so make sure that the MC sample you are using is the + // right one (check all these quantities in your files and compare with + // is written below. In some theta bins, there are two different + // productions, with different energy limits but with the same spectral + // slope. We account for this when calculating the original number of + // events in each energy bin. + // + + TH1D &histsel = *fHistSel->ProjectionX(); + + TH1D histall; + + TAxis &xaxis = *histsel.GetXaxis(); + + MH::SetBinning(&histall, &xaxis); + MH::SetBinning(fHistCol, &xaxis); + + Float_t emin1, emax1, emin2, emax2; + Float_t index, expo, k1, k2; + Float_t numevts1, numevts2; + Float_t r1, r2; // Impact parameter range (on ground). + + emin1 = 0; emax1 = 0; emin2 = 0; emax2 = 0; + expo = 0.; k1 = 0.; k2 = 0.; r1 = 0.; r2 = 0.; + numevts1 = 0; numevts2 = 0; + + if (theta > 0.26 && theta < 0.27) // 15 degrees + { + r1 = 0.; + r2 = 250.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emin2 = 400.; + emax2 = 30000.; + numevts1 = 4000.; + numevts2 = 25740.; + } + else if (theta > 0.34 && theta < 0.35) // 20 degrees + { + r1 = 0.; + r2 = 263.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emin2 = 400.; + emax2 = 30000.; + numevts1 = 6611.; + numevts2 = 24448.; + } + else if (theta > 0.47 && theta < 0.48) // 27 degrees + { + r1 = 0.; + r2 = 290.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emin2 = 400.; + emax2 = 30000.; + numevts1 = 4000.; + numevts2 = 26316.; + } + else if (theta > 0.55 && theta < 0.56) // 32 degrees + { + r1 = 0.; + r2 = 350.; //meters + emin1 = 300.; + emax1 = 30000.; // Energies in GeV. + numevts1 = 33646.; + } + else if (theta > 0.68 && theta < 0.69) // 39 degrees + { + r1 = 0.; + r2 = 380.; //meters + emin1 = 300.; + emax1 = 30000.; // Energies in GeV. + numevts1 = 38415.; + } + else if (theta > 0.78 && theta < 0.79) // 45 degrees + { + r1 = 0.; + r2 = 565.; //meters + emin1 = 300.; + emax1 = 50000.; // Energies in GeV. + numevts1 = 30197.; + } + + index = 1.5; // Differential spectral Index. + expo = 1.-index; + k1 = numevts1 / (pow(emax1,expo) - pow(emin1,expo)); + k2 = numevts2 / (pow(emax2,expo) - pow(emin2,expo)); + + const Int_t nbinx = xaxis.GetNbins(); + + for (Int_t i=1; i<=nbinx; i++) + { + const Float_t e1 = histall.GetBinLowEdge(i); + const Float_t e2 = histall.GetBinLowEdge(i+1); + + if (e1 < emin1 || e2 > emax2) + continue; + + Float_t events; + + if (e2emin2) + events = k2 * (pow(e2, expo) - pow(e1, expo)); + else + events = + k1 * (pow(emax1, expo) - pow(e1, expo))+ + k2 * (pow(e2, expo) - pow(emin2, expo)); + + histall.SetBinContent(i, events); + histall.SetBinError(i, sqrt(events)); + } + + // ----------------------------------------------------------- + + const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1); + + for (Int_t ix=1; ix<=nbinx; ix++) + { + 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 err = sqrt((1.-eff)*Ns)/Na; + + const Float_t area = dr * cos(theta); + + fHistCol->SetBinContent(ix, eff*area); + fHistCol->SetBinError(ix, err*area); + } + + delete &histsel; + + SetReadyToSave(); +} Index: trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h (revision 1821) +++ trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h (revision 1821) @@ -0,0 +1,39 @@ +#ifndef MARS_MHMcCT1CollectionArea +#define MARS_MHMcCT1CollectionArea + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class TH1D; +class TH2D; + + +class MHMcCT1CollectionArea : public MParContainer +{ +private: + TH2D *fHistAll; //! all simulated showers + TH2D *fHistSel; //! the selected showers + + TH1D *fHistCol; // the collection area + +public: + MHMcCT1CollectionArea(const char *name=NULL, const char *title=NULL); + ~MHMcCT1CollectionArea(); + + void FillSel(Double_t energy, Double_t radius); + + void DrawAll(Option_t *option=""); + void DrawSel(Option_t *option=""); + + const TH1D *GetHist() const { return fHistCol; } + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option="") const; + + void CalcEfficiency(Float_t theta); + + ClassDef(MHMcCT1CollectionArea, 1) // Data Container to calculate Collection Area +}; + +#endif Index: trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhist/Makefile (revision 1820) +++ trunk/MagicSoft/Mars/mhist/Makefile (revision 1821) @@ -67,5 +67,6 @@ MHSigmaPixel.cc \ MHSigmabarTheta.cc \ - MHSigmaTheta.cc + MHSigmaTheta.cc \ + MHMcCT1CollectionArea.cc SRCS = $(SRCFILES) Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.cc (revision 1821) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.cc (revision 1821) @@ -0,0 +1,114 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHMcCT1CollectionAreaCalc +// +// +////////////////////////////////////////////////////////////////////////////// + +#include "MMcCT1CollectionAreaCalc.h" + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MMcEvt.hxx" + +#include "MHMcCT1CollectionArea.h" + +ClassImp(MMcCT1CollectionAreaCalc); + +MMcCT1CollectionAreaCalc::MMcCT1CollectionAreaCalc(const char *input, + const char *name, const char *title) +{ + fName = name ? name : "MMcCT1CollectionAreaCalc"; + fTitle = title ? title : "Task to calculate the collection area"; + + AddToBranchList("MMcEvt.fEnergy"); + AddToBranchList("MMcEvt.fImpact"); + AddToBranchList("MMcEvt.fTelescopeTheta"); +} + +Bool_t MMcCT1CollectionAreaCalc::PreProcess (MParList *pList) +{ + // connect the raw data with this task + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... exit." << endl; + return kFALSE; + } + + fCollArea = (MHMcCT1CollectionArea*)pList->FindCreateObj("MHMcCT1CollectionArea"); + if (!fCollArea) + return kFALSE; + + fTotalNumSimulatedShowers = 0; + fAllEvtsTriggered = kTRUE; + fAllEvtsAtSameTheta = kTRUE; + fTelescopeTheta = -1.; + + return kTRUE; +} + +Bool_t MMcCT1CollectionAreaCalc::Process() +{ + const Double_t energy = fMcEvt->GetEnergy(); + const Double_t impact = fMcEvt->GetImpact()/100.; + + if (fTelescopeTheta < 0.) + fTelescopeTheta = fMcEvt->GetTelescopeTheta(); + + if (fTelescopeTheta != fMcEvt->GetTelescopeTheta()) + { + fAllEvtsAtSameTheta = kFALSE; + } + + fCollArea->FillSel(energy, impact); + + return kTRUE; +} + +Bool_t MMcCT1CollectionAreaCalc::PostProcess() +{ + if ( ! fAllEvtsAtSameTheta ) + { + *fLog << err << dbginf << "ERROR: input data contain events at different zenith angles (not supported)... exiting." << endl; + return kFALSE; + } + // + // do the calculation of the effective area + // + *fLog << inf << "Calculation Collection Area..." << endl; + + fCollArea->CalcEfficiency(fTelescopeTheta); + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.h (revision 1821) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.h (revision 1821) @@ -0,0 +1,37 @@ +#ifndef MARS_MMcCT1CollectionAreaCalc +#define MARS_MMcCT1CollectionAreaCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MParList; +class MMcEvt; +class MHMcCT1CollectionArea; + +class MMcCT1CollectionAreaCalc : public MTask +{ +private: + const MMcEvt *fMcEvt; + + MHMcCT1CollectionArea *fCollArea; + + TString fObjName; + + UInt_t fTotalNumSimulatedShowers; + Bool_t fAllEvtsTriggered; + Bool_t fAllEvtsAtSameTheta; + Float_t fTelescopeTheta; + +public: + MMcCT1CollectionAreaCalc(const char *input=NULL, + const char *name=NULL, const char *title=NULL); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + ClassDef(MMcCT1CollectionAreaCalc, 0) // Task to calculate the collection area histogram +}; + +#endif Index: trunk/MagicSoft/Mars/mmontecarlo/Makefile =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/Makefile (revision 1820) +++ trunk/MagicSoft/Mars/mmontecarlo/Makefile (revision 1821) @@ -31,5 +31,6 @@ MMcThresholdCalc.cc \ MMcTimeGenerate.cc \ - MMcTriggerRateCalc.cc + MMcTriggerRateCalc.cc \ + MMcCT1CollectionAreaCalc.cc SRCS = $(SRCFILES) Index: trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h (revision 1820) +++ trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h (revision 1821) @@ -9,4 +9,5 @@ #pragma link C++ class MMcCollectionAreaCalc+; #pragma link C++ class MMcTriggerRateCalc+; +#pragma link C++ class MMcCT1CollectionAreaCalc+; #endif