Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 1973) +++ /trunk/MagicSoft/Mars/Changelog (revision 1974) @@ -1,10 +1,3 @@ -*-*- END OF LINE -*-*- - - 2003/04/22: Abelardo Moralejo - - * mhist/MHMcEnergyMigration.h - Changed ClassDef(MHMcEnergyMigration, 0) - to ClassDef(MHMcEnergyMigration, 1) to allow writing the - histograms to disk. 2003/04/22: Thomas Bretz @@ -15,4 +8,24 @@ * manalysis/MApplyPadding.[h,cc]: - removed + + * mhist/MHMcCT1CollectionArea.[h,cc], + mhist/MHMcCollectionArea.[h,cc], mhist/MHMcDifRate.[h,cc], + mhist/MHMcEfficiency.[h,cc], mhist/MHMcEfficiencyEnergy.[h,cc], + mhist/MHMcEnergy.[h,cc], mhist/MHMcEnergyImpact.[h,cc], + mhist/MHMcEnergyMigration.[h,cc], mhist/MHMcIntRate.[h,cc], + mhist/MHMcRate.[h,cc], mhist/MHMcTriggerLvl2.[h,cc]: + - moved to new directory mhistmc + + * mhistmc, mhistmc/Makefile, mhistmc/HistMcLinkDef.h: + - new + + * mhist/Makefile, mhist/HistLinkDef.h: + - removed MHMc* + + * mhistmc/MHMcCollectionArea.[h,cc]: + - changed to fit requirements of MStatusDisplay + + * Makefile: + - added mhistmc @@ -36,4 +49,9 @@ - Removed. It was not necessary, since the same job can be done with the already existing task MFillH. + + * mhist/MHMcEnergyMigration.h + - Changed ClassDef(MHMcEnergyMigration, 0) + to ClassDef(MHMcEnergyMigration, 1) to allow writing the + histograms to disk. Index: /trunk/MagicSoft/Mars/Makefile =================================================================== --- /trunk/MagicSoft/Mars/Makefile (revision 1973) +++ /trunk/MagicSoft/Mars/Makefile (revision 1974) @@ -46,4 +46,5 @@ mmontecarlo \ mhist \ + mhistmc \ mdata \ mtools \ Index: /trunk/MagicSoft/Mars/mhist/HistLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 1973) +++ /trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 1974) @@ -45,17 +45,3 @@ #pragma link C++ class MHRanForestGini+; -#pragma link C++ class MHMcRate+; -#pragma link C++ class MHMcDifRate+; -#pragma link C++ class MHMcIntRate+; -#pragma link C++ class MHMcEfficiency+; -#pragma link C++ class MHMcEfficiencyImpact+; -#pragma link C++ class MHMcEfficiencyEnergy+; -#pragma link C++ class MHMcEnergy+; -#pragma link C++ class MHMcEnergyImpact+; -#pragma link C++ class MHMcCollectionArea+; -#pragma link C++ class MHMcEnergyMigration+; -#pragma link C++ class MHMcCT1CollectionArea+; - -#pragma link C++ class MHMcTriggerLvl2+; - #endif Index: /trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mhist/Makefile (revision 1973) +++ /trunk/MagicSoft/Mars/mhist/Makefile (revision 1974) @@ -22,5 +22,5 @@ # connect the include files defined in the config.mk file # -INCLUDES = -I. -I../mbase -I../mraw -I../manalysis -I../mmc \ +INCLUDES = -I. -I../mbase -I../mraw -I../manalysis \ -I../mgui -I../mgeom -I../mdata -I../mfilter -I../mimage \ -I../mmain @@ -43,5 +43,4 @@ MHEnergyTime.cc \ MHEnergyTheta.cc \ - MHMcCollectionArea.cc \ MHAlphaEnergyTime.cc \ MHAlphaEnergyTheta.cc \ @@ -54,22 +53,11 @@ MHRanForest.cc \ MHRanForestGini.cc \ - MHMcEnergy.cc \ - MHMcEfficiency.cc \ - MHMcEfficiencyImpact.cc \ - MHMcEfficiencyEnergy.cc \ - MHMcEnergyImpact.cc \ - MHMcRate.cc \ - MHMcIntRate.cc \ - MHMcDifRate.cc \ MHThetabarTime.cc \ MHThetabarTheta.cc \ - MHMcEnergyMigration.cc \ MHGamma.cc \ MHFlux.cc \ MHSigmaPixel.cc \ MHSigmabarTheta.cc \ - MHSigmaTheta.cc \ - MHMcCT1CollectionArea.cc \ - MHMcTriggerLvl2.cc + MHSigmaTheta.cc SRCS = $(SRCFILES) Index: /trunk/MagicSoft/Mars/mhistmc/HistMcIncl.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/HistMcIncl.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/HistMcIncl.h (revision 1974) @@ -0,0 +1,6 @@ +#ifndef __CINT__ + +#include +#include + +#endif // __CINT__ Index: /trunk/MagicSoft/Mars/mhistmc/HistMcLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/HistMcLinkDef.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/HistMcLinkDef.h (revision 1974) @@ -0,0 +1,20 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MHMcRate+; +#pragma link C++ class MHMcDifRate+; +#pragma link C++ class MHMcIntRate+; +#pragma link C++ class MHMcEfficiency+; +#pragma link C++ class MHMcEfficiencyImpact+; +#pragma link C++ class MHMcEfficiencyEnergy+; +#pragma link C++ class MHMcEnergy+; +#pragma link C++ class MHMcEnergyImpact+; +#pragma link C++ class MHMcCollectionArea+; +#pragma link C++ class MHMcEnergyMigration+; +#pragma link C++ class MHMcCT1CollectionArea+; +#pragma link C++ class MHMcTriggerLvl2+; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc (revision 1974) @@ -0,0 +1,368 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// 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, Int_t nbins, Axis_t minEnergy, Axis_t maxEnergy) +{ + // + // nbins, minEnergy, maxEnergy defaults: + // we set the energy range from 100 Gev to 30000 GeV (in log, 3.5 orders + // of magnitude) and for each order we take 10 subdivisions --> 35 xbins + // we set the theta range from 12.5 to 48 deg, with 6 bins (the latter + // choice has been done to make the bin centers as close as possible to + // the actual zenith angles in the CT1 MC sample). + // + + fName = name ? name : "MHMcCT1CollectionArea"; + fTitle = title ? title : "Collection Area vs. log10 Energy"; + + fHistAll = new TH2D; + fHistSel = new TH2D; + fHistCol = new TH2D; + + SetBins(nbins, minEnergy, maxEnergy); + + fHistCol->SetName(fName); + fHistAll->SetName("AllEvents"); + fHistSel->SetName("SelectedEvents"); + + fHistCol->SetTitle(fTitle); + fHistAll->SetTitle("All showers - Theta vs log10 Energy distribution"); + fHistSel->SetTitle("Selected showers - Theta vs log10 Energy distribution"); + + fHistAll->SetDirectory(NULL); + fHistSel->SetDirectory(NULL); + fHistCol->SetDirectory(NULL); + + fHistAll->SetXTitle("log10 E [GeV]"); + fHistAll->SetYTitle("theta [deg]"); + fHistAll->SetZTitle("N"); + + fHistSel->SetXTitle("log10 E [GeV]"); + fHistSel->SetYTitle("theta [deg]"); + fHistSel->SetZTitle("N"); + + fHistCol->SetXTitle("log10 E [GeV]"); + fHistCol->SetYTitle("theta [deg]"); + fHistCol->SetZTitle("A [m^{2}]"); +} + +// -------------------------------------------------------------------------- +// +// Set binning of histograms. Binning of energy axis can be changed, that +// of theta axis is always the same (to match the CT1 MC sample theta dist): +// +void MHMcCT1CollectionArea::SetBins(Int_t nbins, Axis_t minEnergy, Axis_t maxEnergy) +{ + MBinning binsx; + binsx.SetEdges(nbins, minEnergy, maxEnergy); + + MBinning binsy; + const Double_t yedge[7] = {12.5, 17.5, 23.5, 29.5, 35.5, 42., 48.}; + const TArrayD yed(7,yedge); + binsy.SetEdges(yed); + + MH::SetBinning(fHistAll, &binsx, &binsy); + MH::SetBinning(fHistSel, &binsx, &binsy); + MH::SetBinning(fHistCol, &binsx, &binsy); +} + + +// -------------------------------------------------------------------------- +// +// 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 theta) +{ + fHistSel->Fill(log10(energy), theta); +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram with all showers +// +void MHMcCT1CollectionArea::DrawAll(Option_t* option) +{ + if (!gPad) + MH::MakeDefCanvas(fHistAll); + + fHistAll->Draw(option); + + 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->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 behaviour of DrawClone + // + gROOT->SetSelectedPad(NULL); + + fHistCol->DrawCopy(option); + + c->Modified(); + c->Update(); + + return c; +} + +void MHMcCT1CollectionArea::Draw(Option_t* option) +{ + if (!gPad) + MH::MakeDefCanvas(fHistCol); + + fHistCol->Draw(option); + + gPad->Modified(); + gPad->Update(); +} + +// +// Calculate the Efficiency (collection area) for the CT1 MC sample +// and set the 'ReadyToSave' flag +// +void MHMcCT1CollectionArea::CalcEfficiency() +{ + // + // 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. + // + + for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++) + { + // This theta is not exactly the one of the MC events, just about + // the same: + Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin); + + 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 > 14 && theta < 16) // 15 deg + { + r1 = 0.; + r2 = 250.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emin2 = 400.; + emax2 = 30000.; + numevts1 = 4000.; + numevts2 = 25740.; + } + else if (theta > 20 && theta < 21) // 20.5 deg + { + r1 = 0.; + r2 = 263.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emin2 = 400.; + emax2 = 30000.; + numevts1 = 6611.; + numevts2 = 24448.; + } + else if (theta > 26 && theta < 27) // 26.5 degrees + { + r1 = 0.; + r2 = 290.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emax2 = emax1; emin2 = 400.; + emax2 = 30000.; + numevts1 = 4000.; + numevts2 = 26316.; + } + else if (theta > 32 && theta < 33) // 32.5 degrees + { + r1 = 0.; + r2 = 350.; //meters + emin1 = 300.; + emax1 = 30000.; // Energies in GeV. + emax2 = emax1; + numevts1 = 33646.; + } + else if (theta > 38 && theta < 39) // 38.75 degrees + { + r1 = 0.; + r2 = 380.; //meters + emin1 = 300.; + emax1 = 30000.; // Energies in GeV. + emax2 = emax1; + numevts1 = 38415.; + } + else if (theta > 44 && theta < 46) // 45 degrees + { + r1 = 0.; + r2 = 565.; //meters + emin1 = 300.; + emax1 = 50000.; // Energies in GeV. + emax2 = emax1; + 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)); + + for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++) + { + const Float_t e1 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i)); + const Float_t e2 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i+1)); + + if (e1 < emin1 || e2 > emax2) + continue; + + Float_t events; + + if (e2 <= emax1) + events = k1 * (pow(e2, expo) - pow(e1, expo)); + else if (e1 >= emin2) + events = k2 * (pow(e2, expo) - pow(e1, expo)); + else + events = + k1 * (pow(emax1, expo) - pow(e1, expo))+ + k2 * (pow(e2, expo) - pow(emin2, expo)); + + fHistAll->SetBinContent(i, thetabin, events); + fHistAll->SetBinError(i, thetabin, sqrt(events)); + } + + // ----------------------------------------------------------- + + const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1); + + for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++) + { + const Float_t Na = fHistAll->GetBinContent(ix,thetabin); + + if (Na <= 0) + { + // + // If energy is large, this case means that no or very few events + // were generated at this energy bin. In this case we assign it + // the effective area of the bin below it in energy. If energy is + // below 1E4, it means that no events triggered -> eff area = 0 + // + + if (fHistSel->GetXaxis()->GetBinLowEdge(ix) > 4.) + { + fHistCol->SetBinContent(ix, thetabin, fHistCol->GetBinContent(ix-1, thetabin)); + fHistCol->SetBinError(ix, thetabin, fHistCol->GetBinError(ix-1, thetabin)); + } + continue; + } + + const Float_t Ns = fHistSel->GetBinContent(ix,thetabin); + + // 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*TMath::Pi()/180.); + + fHistCol->SetBinContent(ix, thetabin, eff*area); + fHistCol->SetBinError(ix, thetabin, err*area); + + } + } + + SetReadyToSave(); +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.h (revision 1974) @@ -0,0 +1,44 @@ +#ifndef MARS_MHMcCT1CollectionArea +#define MARS_MHMcCT1CollectionArea + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class TH2D; + + +class MHMcCT1CollectionArea : public MParContainer +{ +private: + TH2D *fHistAll; // all simulated showers + TH2D *fHistSel; // the selected showers + TH2D *fHistCol; // the collection area + void SetBins(Int_t nbins, Axis_t minEnergy, Axis_t maxEnergy); + +public: + MHMcCT1CollectionArea(const char *name=NULL, const char *title=NULL, Int_t nbins=30, Axis_t minEnergy=2., Axis_t maxEnergy=5.); + ~MHMcCT1CollectionArea(); + + void FillSel(Double_t energy, Double_t radius); + + void DrawAll(Option_t *option=""); + void DrawSel(Option_t *option=""); + + const TH2D *GetHist() const { return fHistCol; } + const TH2D *GetHAll() const { return fHistAll; } + const TH2D *GetHSel() const { return fHistSel; } + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option="") const; + + void CalcEfficiency(); + + ClassDef(MHMcCT1CollectionArea, 1) // Data Container to calculate Collection Area +}; + +#endif + + + + Index: /trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc (revision 1974) @@ -0,0 +1,382 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHMcCollectionArea // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHMcCollectionArea.h" + +#include +#include + +#include "MH.h" +#include "MBinning.h" +#include "MHMcEfficiency.h" +#include "MHMcEnergyImpact.h" + +ClassImp(MHMcCollectionArea); + +// -------------------------------------------------------------------------- +// +// Creates the three necessary histograms: +// - selected showers (input) +// - all showers (input) +// - collection area (result) +// +MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title) +{ + // initialize the histogram for the distribution r vs E + // + // we set the energy range from 10 Gev to 20000 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 : "Collection Area vs. Energy"; + + MBinning binsx; + MBinning binsy; + binsx.SetEdgesLog(50, 10, 20000); + 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 +// +MHMcCollectionArea::~MHMcCollectionArea() +{ + delete fHistAll; + delete fHistSel; + delete fHistCol; +} + +// -------------------------------------------------------------------------- +// +// Fill data into the histogram which contains all showers +// +void MHMcCollectionArea::FillAll(Double_t energy, Double_t radius) +{ + fHistAll->Fill(energy, radius); +} + +// -------------------------------------------------------------------------- +// +// 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(); +} + +// -------------------------------------------------------------------------- +// +// 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 *MHMcCollectionArea::DrawClone(Option_t* opt) const +{ + return MH::DrawClone(opt); +/* + 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 MHMcCollectionArea::Draw(Option_t* option) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + // if (!gPad) + // MH::MakeDefCanvas(fHistCol); + + 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(); + + const Double_t max = psel.GetYaxis()->GetXmax(); + + fHistCol->Divide(&psel, &pall, TMath::Pi()*max*max, 1); + fHistCol->SetEntries(hsel.GetEntries()); + + delete &psel; + delete &pall; + + SetReadyToSave(); +} + +// -------------------------------------------------------------------------- +// +// Calculate the Efficiency (collection area) and set the 'ReadyToSave' +// flag +// +void MHMcCollectionArea::CalcEfficiency(UInt_t numevts, Float_t angle) +{ + // Here we estimate the total number of showers in each energy bin + // known the energy range and spectral index of the generated showers + // (set now by hand since the info is not available in the header!) + + TH1D &histsel = *fHistSel->ProjectionX(); + + TH1D histall; + + TAxis &xaxis = *histsel.GetXaxis(); + + MH::SetBinning(&histall, &xaxis); + MH::SetBinning(fHistCol, &xaxis); + + // TO BE FIXED! In forthcoming camera versions emin, emax and index + // will be available in a run header, and should be read from it. + + const Float_t emin = 10.; + const Float_t emax = 30000.; // Energies in GeV. + const Float_t index = 2.6; // Differential spectral Index. + + const Float_t expo = 1.-index; + + const Float_t k = (Float_t)numevts / (pow(emax,expo) - pow(emin,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 < emin || e2 > emax) + continue; + + const Float_t events = k * (pow(e2, expo) - pow(e1, expo)); + + histall.SetBinContent(i, events); + histall.SetBinError(i, sqrt(events)); + + } + + // ----------------------------------------------------------- + + // Impact parameter range: + const Float_t r1 = 0; + const Float_t r2 = 400; + + const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1); + + angle *= TMath::Pi()/180; + + 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(angle); + + fHistCol->SetBinContent(ix, eff*area); + fHistCol->SetBinError(ix, err*area); + } + + 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 err = h.GetCellError(ix, iy); + + colA += eff*A; + errA += A*A * err*err; + } + + fHistCol->SetBinContent(ix, colA); + fHistCol->SetBinError(ix, sqrt(errA)); + } + + SetReadyToSave(); +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h (revision 1974) @@ -0,0 +1,49 @@ +#ifndef MARS_MHMcCollectionArea +#define MARS_MHMcCollectionArea + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TH2D; + +class MHMcEfficiency; +class MHMcEnergyImpact; + +class MHMcCollectionArea : public MH +{ +private: + TH2D *fHistAll; //! all simulated showers + TH2D *fHistSel; //! the selected showers + + TH1D *fHistCol; // the collection area + + void Calc(TH2D &hsel, TH2D &hall); + +public: + 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 DrawAll(Option_t *option=""); + void DrawSel(Option_t *option=""); + + const TH1D *GetHist() { return fHistCol; } + const TH1D *GetHist() const { return fHistCol; } + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option="") const; + + void CalcEfficiency(); + void CalcEfficiency(UInt_t allevts, Float_t theta); + + void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall); + void Calc(const MHMcEfficiency &heff); + + ClassDef(MHMcCollectionArea, 1) // Data Container to calculate Collection Area +}; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcDifRate.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcDifRate.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcDifRate.cc (revision 1974) @@ -0,0 +1,183 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Lopez 05/2001 +! Author(s): Thomas Bretz 05/2001 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMcDifRate +// +// This class holds the information (histogram and fit function) +// about the energy threshold for a particular trigger condition. +// +//////////////////////////////////////////////////////////////////////////// +#include "MHMcDifRate.h" + +#include + +#include + +#include "MH.h" +#include "MBinning.h" + +#include "MHMcCollectionArea.h" + +ClassImp(MHMcDifRate); + +// ------------------------------------------------------------------------- +// +// Default Constructor. +// +MHMcDifRate::MHMcDifRate(const char *name, const char *title) + : fHist() +{ + fName = name ? name : "MHMcDifRate"; + fTitle = title ? title : "Differential Trigger Rate"; + + // - we initialize the histogram + // - we have to give diferent names for the diferent histograms because + // root don't allow us to have diferent histograms with the same name + + fHist.SetName(fName); + fHist.SetTitle(fTitle); + + + fHist.SetDirectory(NULL); + + fHist.SetXTitle("E [GeV]"); + fHist.SetYTitle("dR/dE [Hz/GeV]"); +} + +void MHMcDifRate::SetName(const char *name) +{ + fName = name; + fHist.SetName(name); + fHist.SetDirectory(NULL); +} + +void MHMcDifRate::SetTitle(const char *title) +{ + fTitle = title; + fHist.SetTitle(title); +} + +//------------------------------------------------------------------------- +// +// Defualt Destructor +// +MHMcDifRate::~MHMcDifRate() +{ +} + +// ------------------------------------------------------------------------ +// +// Drawing function. It creates its own canvas. +// +void MHMcDifRate::Draw(Option_t *option) +{ + if (!gPad) + MH::MakeDefCanvas(&fHist); + + gPad->SetLogx(); + + fHist.Draw(option); + + gPad->Modified(); + gPad->Update(); +} + +TObject *MHMcDifRate::DrawClone(Option_t *option) const +{ + TCanvas *c = MH::MakeDefCanvas(&fHist); + + c->SetLogx(); + + // + // This is necessary to get the expected bahviour of DrawClone + // + gROOT->SetSelectedPad(NULL); + + ((TH1D&)fHist).DrawCopy(option); + + c->Modified(); + c->Update(); + + return c; +} +/* +void MHMcDifRate::Calc(const TH2D &hsel, const TH2D &hall) +{ + // + // Set the binning from the two axis of one of the two histograms + // + MH::SetBinning(&fHist, ((TH2D&)hsel).GetXaxis(), ((TH2D&)hsel).GetYaxis()); + + // + // This is necessary to initialize thze error calculation correctly + // (Nothing important: The histogram set the size of its internal + // array storing the errors to the correct size) + // + fHist.Sumw2(); + + // + // Calculate the efficiency by dividing the number of selected + // (eg. triggered) showers by the number of all showers per bin. + // Both histograms are weighted with weight 1, and for the error + // calculation we assume a binomial error calculation. + // + fHist.Divide((TH2D*)&hsel, (TH2D*)&hall, 1, 1, "B"); + + SetReadyToSave(); +} +*/ + +// -------------------------------------------------------------------------- +// +// Calculate the DifRate and set the 'ReadyToSave' flag. +// The DifRate is calculated as the number of selected showers +// (eg. triggered ones) divided by the number of all showers. +// For error calculation Binomial errors are computed. +// +void MHMcDifRate::Calc(const MHMcCollectionArea &cola, const TF1 &spect) +{ + /*const*/ TH1D &hcol = (TH1D&)*cola.GetHist(); + + fHist.Reset(); + + // + // Set the binning from the two axis of one of the two histograms + // + MH::SetBinning(&fHist, hcol.GetXaxis()); + + // + // This is necessary to initialize thze error calculation correctly + // (Nothing important: The histogram set the size of its internal + // array storing the errors to the correct size) + // + fHist.Sumw2(); + + fHist.Add(&hcol); + fHist.Multiply((TF1*)&spect); + + SetReadyToSave(); +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcDifRate.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcDifRate.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcDifRate.h (revision 1974) @@ -0,0 +1,44 @@ +#ifndef MARS_MHMcDifRate +#define MARS_MHMcDifRate + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +#ifndef ROOT_TH1 +#include +#endif + +class TF1; +class MMcEvt; +class MHMcCollectionArea; + +class MHMcDifRate : public MParContainer +{ +private: + + TH1D fHist; // histogram with the logarith of the energy + + // const MMcEvt *fMcEvt; //! container to fill histogram from + +public: + MHMcDifRate(const char *name=NULL, const char *title=NULL); + ~MHMcDifRate(); + + void SetName(const char *name); + void SetTitle(const char *title); + + const TH1D *GetHist() { return &fHist; } + const TH1D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t* option = ""); + TObject *DrawClone(Option_t* option = "") const; + +// void Calc(const TH2D &hsel, const TH2D &hall); + void Calc(const MHMcCollectionArea &cola, const TF1 &spect); + + ClassDef(MHMcDifRate, 1) // Histogram container for differential trigger rate +}; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiency.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiency.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiency.cc (revision 1974) @@ -0,0 +1,173 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Lopez 05/2001 +! Author(s): Thomas Bretz 05/2001 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMcEfficiency +// +// This class holds the information (histogram and fit function) +// about the energy threshold for a particular trigger condition. +// +//////////////////////////////////////////////////////////////////////////// +#include "MHMcEfficiency.h" + +#include + +#include +#include + +#include "MH.h" +#include "MBinning.h" + +#include "MHMcEnergyImpact.h" + +ClassImp(MHMcEfficiency); + +// ------------------------------------------------------------------------- +// +// Default Constructor. +// +MHMcEfficiency::MHMcEfficiency(const char *name, const char *title) + : fHist() +{ + fName = name ? name : "MHMcEfficiency"; + fTitle = title ? title : "Trigger Efficieny (Energy-Impact parameter plane)"; + + // - we initialize the histogram + // - we have to give diferent names for the diferent histograms because + // root don't allow us to have diferent histograms with the same name + + fHist.SetName(fName); + fHist.SetTitle(fTitle); + + fHist.SetDirectory(NULL); + + fHist.SetXTitle("E [GeV]"); + fHist.SetYTitle("r [m]"); + fHist.SetZTitle("Trig. Eff. [1]"); + + + MBinning binsx; + binsx.SetEdgesLog(10, 1, 100000); // [GeV] + + MBinning binsy; + binsy.SetEdges(9, 0, 450); // [m] + MH::SetBinning(&fHist, &binsx, &binsy); +} + +void MHMcEfficiency::SetName(const char *name) +{ + fName = name; + fHist.SetName(name); + fHist.SetDirectory(NULL); +} + +void MHMcEfficiency::SetTitle(const char *title) +{ + fTitle = title; + fHist.SetTitle(title); +} + +//------------------------------------------------------------------------- +// +// Defualt Destructor +// +MHMcEfficiency::~MHMcEfficiency() +{ +} + +// ------------------------------------------------------------------------ +// +// Drawing function. It creates its own canvas. +// +void MHMcEfficiency::Draw(Option_t *option) +{ + if (!gPad) + MH::MakeDefCanvas(&fHist); + + gPad->SetLogx(); + + fHist.Draw(option); + + gPad->Modified(); + gPad->Update(); +} + +TObject *MHMcEfficiency::DrawClone(Option_t *option) const +{ + TCanvas *c = MH::MakeDefCanvas(&fHist); + + c->SetLogx(); + + // + // This is necessary to get the expected bahviour of DrawClone + // + gROOT->SetSelectedPad(NULL); + + ((TH2D&)fHist).DrawCopy(option); + + c->Modified(); + c->Update(); + + return c; +} + +void MHMcEfficiency::Calc(const TH2D &hsel, const TH2D &hall) +{ + // + // Set the binning from the two axis of one of the two histograms + // + MH::SetBinning(&fHist, ((TH2D&)hsel).GetXaxis(), ((TH2D&)hsel).GetYaxis()); + + // + // This is necessary to initialize thze error calculation correctly + // (Nothing important: The histogram set the size of its internal + // array storing the errors to the correct size) + // + fHist.Sumw2(); + + // + // Calculate the efficiency by dividing the number of selected + // (eg. triggered) showers by the number of all showers per bin. + // Both histograms are weighted with weight 1, and for the error + // calculation we assume a binomial error calculation. + // + fHist.Divide((TH2D*)&hsel, (TH2D*)&hall, 1, 1, "B"); + + SetReadyToSave(); +} + +// -------------------------------------------------------------------------- +// +// Calculate the Efficiency and set the 'ReadyToSave' flag. +// The Efficiency is calculated as the number of selected showers +// (eg. triggered ones) divided by the number of all showers. +// For error calculation Binomial errors are computed. +// +void MHMcEfficiency::Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall) +{ + Calc(*mcsel.GetHist(), *mcall.GetHist()); + +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiency.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiency.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiency.h (revision 1974) @@ -0,0 +1,43 @@ +#ifndef MARS_MHMcEfficiency +#define MARS_MHMcEfficiency + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +#ifndef ROOT_TH2 +#include +#endif + +class MMcEvt; +class MHMcEnergyImpact; + +class MHMcEfficiency : public MParContainer +{ +private: + + TH2D fHist; // histogram with the logarith of the energy + + // const MMcEvt *fMcEvt; //! container to fill histogram from + +public: + MHMcEfficiency(const char *name=NULL, const char *title=NULL); + ~MHMcEfficiency(); + + void SetName(const char *name); + void SetTitle(const char *title); + + const TH2D *GetHist() { return &fHist; } + const TH2D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t* option = ""); + TObject *DrawClone(Option_t* option = "") const; + + void Calc(const TH2D &hsel, const TH2D &hall); + void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall); + + ClassDef(MHMcEfficiency, 1) // Histogram container for 2D trigger efficiency +}; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyEnergy.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyEnergy.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyEnergy.cc (revision 1974) @@ -0,0 +1,171 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Lopez 05/2001 +! Author(s): Thomas Bretz 05/2001 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMcEfficiencyEnergy +// +// This class holds the information (histogram and fit function) +// about the energy threshold for a particular trigger condition. +// +//////////////////////////////////////////////////////////////////////////// +#include "MHMcEfficiencyEnergy.h" + +#include + +#include +#include + +#include "MH.h" +#include "MBinning.h" + +#include "MHMcEnergyImpact.h" + +ClassImp(MHMcEfficiencyEnergy); + +// ------------------------------------------------------------------------- +// +// Default Constructor. +// +MHMcEfficiencyEnergy::MHMcEfficiencyEnergy(const char *name, const char *title) + : fHist() +{ + fName = name ? name : "MHMcEfficiencyEnergy"; + fTitle = title ? title : "Trigger Efficieny vs. Energy"; + + // - we initialize the histogram + // - we have to give diferent names for the diferent histograms because + // root don't allow us to have diferent histograms with the same name + + fHist.SetName(fName); + fHist.SetTitle(fTitle); + + fHist.SetDirectory(NULL); + + fHist.SetXTitle("E [GeV]"); + fHist.SetYTitle("Trig. Eff. [1]"); + + MBinning binse; + binse.SetEdgesLog(10, 1, 100000); // [GeV] + MH::SetBinning(&fHist, &binse); +} + +void MHMcEfficiencyEnergy::SetName(const char *name) +{ + fName = name; + fHist.SetName(name); + fHist.SetDirectory(NULL); +} + +void MHMcEfficiencyEnergy::SetTitle(const char *title) +{ + fTitle = title; + fHist.SetTitle(title); +} + +//------------------------------------------------------------------------- +// +// Defualt Destructor +// +MHMcEfficiencyEnergy::~MHMcEfficiencyEnergy() +{ +} + +// ------------------------------------------------------------------------ +// +// Drawing function. It creates its own canvas. +// +void MHMcEfficiencyEnergy::Draw(Option_t *option) +{ + if (!gPad) + MH::MakeDefCanvas(&fHist); + + gPad->SetLogx(); + + fHist.Draw(option); + + gPad->Modified(); + gPad->Update(); +} + +TObject *MHMcEfficiencyEnergy::DrawClone(Option_t *option) const +{ + TCanvas *c = MH::MakeDefCanvas(&fHist); + + c->SetLogx(); + + // + // This is necessary to get the expected bahviour of DrawClone + // + gROOT->SetSelectedPad(NULL); + + ((TH1D&)fHist).DrawCopy(option); + + c->Modified(); + c->Update(); + + return c; +} + +void MHMcEfficiencyEnergy::Calc(const TH2D &hsel, const TH2D &hall) +{ + // + // Set the binning from the two axis of one of the two histograms + // + MH::SetBinning(&fHist, ((TH2D&)hsel).GetXaxis()); + + // + // This is necessary to initialize thze error calculation correctly + // (Nothing important: The histogram set the size of its internal + // array storing the errors to the correct size) + // + fHist.Sumw2(); + + // + // Calculate the efficiency by dividing the number of selected + // (eg. triggered) showers by the number of all showers per bin. + // Both histograms are weighted with weight 1, and for the error + // calculation we assume a binomial error calculation. + // + TH1D &tsel = *((TH2D&)hsel).ProjectionX(); + TH1D &tall = *((TH2D&)hall).ProjectionX(); + fHist.Divide(&tsel, &tall, 1, 1); + delete &tsel; + delete &tall; + + SetReadyToSave(); +} + +// -------------------------------------------------------------------------- +// +// Calculate the EfficiencyEnergy and set the 'ReadyToSave' flag. +// The EfficiencyEnergy is calculated as the number of selected showers +// (eg. triggered ones) divided by the number of all showers. +// For error calculation Binomial errors are computed. +// +void MHMcEfficiencyEnergy::Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall) +{ + Calc(*mcsel.GetHist(), *mcall.GetHist()); +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyEnergy.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyEnergy.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyEnergy.h (revision 1974) @@ -0,0 +1,43 @@ +#ifndef MARS_MHMcEfficiencyEnergy +#define MARS_MHMcEfficiencyEnergy + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +#ifndef ROOT_TH2 +#include +#endif + +class MMcEvt; +class MHMcEnergyImpact; + +class MHMcEfficiencyEnergy : public MParContainer +{ +private: + + TH1D fHist; // histogram with the logarith of the energy + + // const MMcEvt *fMcEvt; //! container to fill histogram from + +public: + MHMcEfficiencyEnergy(const char *name=NULL, const char *title=NULL); + ~MHMcEfficiencyEnergy(); + + void SetName(const char *name); + void SetTitle(const char *title); + + const TH1D *GetHist() { return &fHist; } + const TH1D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t* option = ""); + TObject *DrawClone(Option_t* option = "") const; + + void Calc(const TH2D &hsel, const TH2D &hall); + void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall); + + ClassDef(MHMcEfficiencyEnergy, 1) // Histogram container for 1D trigger efficiency in energy +}; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyImpact.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyImpact.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyImpact.cc (revision 1974) @@ -0,0 +1,171 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Lopez 05/2001 +! Author(s): Thomas Bretz 05/2001 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMcEfficiencyImpact +// +// This class holds the information (histogram and fit function) +// about the energy threshold for a particular trigger condition. +// +//////////////////////////////////////////////////////////////////////////// +#include "MHMcEfficiencyImpact.h" + +#include + +#include + +#include "MH.h" +#include "MBinning.h" + +#include "MHMcEnergyImpact.h" + +ClassImp(MHMcEfficiencyImpact); + +// ------------------------------------------------------------------------- +// +// Default Constructor. +// +MHMcEfficiencyImpact::MHMcEfficiencyImpact(const char *name, const char *title) + : fHist() +{ + fName = name ? name : "MHMcEfficiencyImpact"; + fTitle = title ? title : "Trigger Efficieny vs. Impact"; + + // - we initialize the histogram + // - we have to give diferent names for the diferent histograms because + // root don't allow us to have diferent histograms with the same name + + fHist.SetName(fName); + fHist.SetTitle(fTitle); + + fHist.SetDirectory(NULL); + + fHist.SetXTitle("r [m]"); + fHist.SetYTitle("Trig. Eff. [1]"); + + + MBinning binsr; + binsr.SetEdges(9, 0, 450); // [m] + MH::SetBinning(&fHist, &binsr); +} + +void MHMcEfficiencyImpact::SetName(const char *name) +{ + fName = name; + fHist.SetName(name); + fHist.SetDirectory(NULL); +} + +void MHMcEfficiencyImpact::SetTitle(const char *title) +{ + fTitle = title; + fHist.SetTitle(title); +} + +//------------------------------------------------------------------------- +// +// Defualt Destructor +// +MHMcEfficiencyImpact::~MHMcEfficiencyImpact() +{ +} + +// ------------------------------------------------------------------------ +// +// Drawing function. It creates its own canvas. +// +void MHMcEfficiencyImpact::Draw(Option_t *option) +{ + if (!gPad) + MH::MakeDefCanvas(&fHist); + + gPad->SetLogx(); + + fHist.Draw(option); + + gPad->Modified(); + gPad->Update(); +} + +TObject *MHMcEfficiencyImpact::DrawClone(Option_t *option) const +{ + TCanvas *c = MH::MakeDefCanvas(&fHist); + + c->SetLogx(); + + // + // This is necessary to get the expected bahviour of DrawClone + // + gROOT->SetSelectedPad(NULL); + + ((TH1D&)fHist).DrawCopy(option); + + c->Modified(); + c->Update(); + + return c; +} + +void MHMcEfficiencyImpact::Calc(const TH2D &hsel, const TH2D &hall) +{ + // + // Set the binning from the two axis of one of the two histograms + // + MH::SetBinning(&fHist, ((TH2D&)hsel).GetYaxis()); + + // + // This is necessary to initialize thze error calculation correctly + // (Nothing important: The histogram set the size of its internal + // array storing the errors to the correct size) + // + fHist.Sumw2(); + + // + // Calculate the efficiency by dividing the number of selected + // (eg. triggered) showers by the number of all showers per bin. + // Both histograms are weighted with weight 1, and for the error + // calculation we assume a binomial error calculation. + // + TH1D &tsel = *((TH2D&)hsel).ProjectionY(); + TH1D &tall = *((TH2D&)hall).ProjectionY(); + fHist.Divide(&tsel, &tall, 1, 1); + delete &tsel; + delete &tall; + + SetReadyToSave(); +} + +// -------------------------------------------------------------------------- +// +// Calculate the EfficiencyImpact and set the 'ReadyToSave' flag. +// The EfficiencyImpact is calculated as the number of selected showers +// (eg. triggered ones) divided by the number of all showers. +// For error calculation Binomial errors are computed. +// +void MHMcEfficiencyImpact::Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall) +{ + Calc(*mcsel.GetHist(), *mcall.GetHist()); +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyImpact.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyImpact.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEfficiencyImpact.h (revision 1974) @@ -0,0 +1,43 @@ +#ifndef MARS_MHMcEfficiencyImpact +#define MARS_MHMcEfficiencyImpact + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +#ifndef ROOT_TH2 +#include +#endif + +class MMcEvt; +class MHMcEnergyImpact; + +class MHMcEfficiencyImpact : public MParContainer +{ +private: + + TH1D fHist; // histogram with the logarith of the energy + + // const MMcEvt *fMcEvt; //! container to fill histogram from + +public: + MHMcEfficiencyImpact(const char *name=NULL, const char *title=NULL); + ~MHMcEfficiencyImpact(); + + void SetName(const char *name); + void SetTitle(const char *title); + + const TH1D *GetHist() { return &fHist; } + const TH1D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t* option = ""); + TObject *DrawClone(Option_t* option = "") const; + + void Calc(const TH2D &hsel, const TH2D &hall); + void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall); + + ClassDef(MHMcEfficiencyImpact, 1) // Histogram container for 1D trigger efficiency in impact +}; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEnergy.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEnergy.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEnergy.cc (revision 1974) @@ -0,0 +1,265 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Lopez 05/2001 +! Author(s): Thomas Bretz 05/2001 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMcEnergy +// +// This class holds the information (histogram and fit function) +// about the energy threshold for a particular trigger condition. +// +//////////////////////////////////////////////////////////////////////////// +#include "MHMcEnergy.h" + +#include +#include + +#include +#include +#include +#include + +#include "MH.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHMcEnergy); + +// ------------------------------------------------------------------------- +// +// Default Constructor. +// +MHMcEnergy::MHMcEnergy(const char *name, const char *title) +{ + fTitle = title ? title : "Container for an energy distribution histogram"; + + // - we initialize the histogram + // - we have to give diferent names for the diferent histograms because + // root don't allow us to have diferent histograms with the same name + + fHist = new TH1F("", "", 20, 0.5, 4.5); + + fHist->SetDirectory(NULL); + fHist->SetXTitle("log(E/GeV)"); + fHist->SetYTitle("dN/dE"); + + SetName(name ? name : "MHMcEnergy"); +} + +// ------------------------------------------------------------------------- +// +// This doesn't only set the name. It tries to get the number from the +// name and creates also name and title of the histogram. +// +// This is necessary for example if a list of such MHMcEnergy histograms +// is created (s. MParList::CreateObjList) +// +void MHMcEnergy::SetName(const char *name) +{ + TString cname(name); + const char *semicolon = strrchr(cname, ';'); + + UInt_t idx = semicolon ? atoi(semicolon+1) : 0; + + fName = cname; + + char text[256]; + if (idx>0) + sprintf(text, "Energy Distribution for trigger condition #%i", idx); + else + sprintf(text, "Energy Distribution"); + + char aux[256]; + strcpy(aux, "Threshold"); + + if (idx>0) + sprintf(aux+strlen(aux), " #%i", idx); + + fHist->SetName(aux); + fHist->SetTitle(text); +} + +//------------------------------------------------------------------------- +// +// Defualt Destructor +// +MHMcEnergy::~MHMcEnergy() +{ + delete fHist; +} + +//-------------------------------------------------------------------------- +// +// Fill the histogram with the log10 of the energy for triggered events. +// +void MHMcEnergy::Fill(Float_t log10E, Float_t w) +{ + fHist->Fill(log10E, w); +} + +// ------------------------------------------------------------------------- +// +// Fitting function +// +void MHMcEnergy::Fit(Axis_t xxmin, Axis_t xxmax) +{ + // + // 0: don't draw the function (it is drawn together with the histogram) + // Q: quiet mode + // + fHist->Fit("gaus", "Q0", "", xxmin, xxmax); + + TF1 *result = fHist->GetFunction("gaus"); + + fThreshold = CalcThreshold(result); + fThresholdErr = CalcThresholdErr(result); + fGaussPeak = CalcGaussPeak(result); + fGaussSigma = CalcGaussSigma(result); +} + +// ------------------------------------------------------------------------ +// +// Helper function for Draw() and DrawClone() which adds some useful +// information to the plot. +// +void MHMcEnergy::DrawLegend() const +{ + char text[256]; + + const Float_t min = fHist->GetMinimum(); + const Float_t max = fHist->GetMaximum(); + const Float_t sum = min+max; + + sprintf(text, "Energy Threshold = %4.1f +- %4.1f GeV", + fThreshold, fThresholdErr); + + TPaveLabel* label = new TPaveLabel(2.2, 0.75*sum, 4.4, 0.90*sum, text); + label->SetFillColor(10); + label->SetTextSize(0.3); + label->SetBit(kCanDelete); + label->Draw(); +} + +// ------------------------------------------------------------------------ +// +// Drawing function. It creates its own canvas. +// +void MHMcEnergy::Draw(Option_t *option) +{ + if (!gPad) + MH::MakeDefCanvas(fHist); + + fHist->Draw(option); + + DrawLegend(); + + gPad->Modified(); + gPad->Update(); +} + +TObject *MHMcEnergy::DrawClone(Option_t *option) const +{ + TCanvas *c = MH::MakeDefCanvas(fHist); + + // + // This is necessary to get the expected bahviour of DrawClone + // + gROOT->SetSelectedPad(NULL); + + fHist->DrawClone(option); + + DrawLegend(); + + c->Modified(); + c->Update(); + + return c; +} + +// -------------------------------------------------------------------------- +// +// Set the number of bins in the histogran. +// +void MHMcEnergy::SetNumBins(Int_t nbins) +{ + fHist->SetBins(nbins, 0.5, 4.5); +} +// -------------------------------------------------------------------------- +// +// Write the threshold and its error in the standard output +// +void MHMcEnergy::Print(Option_t*) const +{ + *fLog << all << "Threshold: " << fThreshold << " +- " << fThresholdErr << endl; +} + +// ------------------------------------------------------------------------- +// +// Return the threshold +// +Float_t MHMcEnergy::CalcThreshold(TF1 *gauss) +{ + const Float_t p1 = gauss->GetParameter(1); + + return pow(10, p1); +} + +// ------------------------------------------------------------------------- +// +// Return the error of the threshold. +// +Float_t MHMcEnergy::CalcThresholdErr(TF1 *gauss) +{ + const Float_t lg10 = log(10); + const Float_t p1 = gauss->GetParameter(1); + const Float_t p1err = gauss->GetParError(1); + + // The error has into accuont the error in the fit + return pow(10, p1) * p1err * lg10; +} + +// ------------------------------------------------------------------------- +// +// Return the peak of the fitted gaussan function. +// +Float_t MHMcEnergy::CalcGaussPeak(TF1 *gauss) +{ + return gauss->GetParameter(1); +} + +// ------------------------------------------------------------------------- +// +// Return the sigma of the fitted gaussan function. +// +Float_t MHMcEnergy::CalcGaussSigma(TF1 *gauss) +{ + return gauss->GetParameter(2); +} + +TH1 *MHMcEnergy::GetHistByName(const TString name) +{ + return fHist; +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEnergy.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEnergy.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEnergy.h (revision 1974) @@ -0,0 +1,58 @@ +#ifndef MARS_MHMcEnergy +#define MARS_MHMcEnergy + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class TH1; +class TH1F; +class TF1; + +class MHMcEnergy : public MParContainer +{ +private: + + TH1F *fHist; // histogram with the logarith of the energy + + Float_t fThreshold; + Float_t fThresholdErr; + Float_t fGaussPeak; + Float_t fGaussSigma; + + Float_t CalcThreshold(TF1 *gauss); + Float_t CalcThresholdErr(TF1 *gauss); + + Float_t CalcGaussPeak(TF1 *gauss); + Float_t CalcGaussSigma(TF1 *gauss); + + void DrawLegend() const; + +public: + + MHMcEnergy(const char *name=NULL, const char *title=NULL); + ~MHMcEnergy(); + + void SetName(const char *name); + + Float_t GetThreshold() const { return fThreshold; } + Float_t GetThresholdErr() const { return fThresholdErr; } + + Float_t GetGaussPeak() const { return fGaussPeak; } + Float_t GetGaussSigma() const { return fGaussSigma; }; + + void Fill(Float_t log10E, Float_t w); + void Fit(Axis_t xxmin, Axis_t xxmax); + void SetNumBins(Int_t nbins = 100); + + TH1 *GetHistByName(const TString name); + + void Draw(Option_t* option = ""); + TObject *DrawClone(Option_t* option = "") const; + + void Print(Option_t* option = NULL) const; + + ClassDef(MHMcEnergy, 1) // Histogram container for montecarlo energy threshold +}; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyImpact.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyImpact.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyImpact.cc (revision 1974) @@ -0,0 +1,168 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Lopez 05/2001 +! Author(s): Thomas Bretz 05/2001 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMcEnergyImpact +// +// This class holds the information (histogram and fit function) +// about the energy threshold for a particular trigger condition. +// +//////////////////////////////////////////////////////////////////////////// +#include "MHMcEnergyImpact.h" + +#include +#include + +#include "MParList.h" +#include "MBinning.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MMcEvt.hxx" + +ClassImp(MHMcEnergyImpact); + +// ------------------------------------------------------------------------- +// +// Default Constructor. +// +MHMcEnergyImpact::MHMcEnergyImpact(const char *name, const char *title) + : fHist() +{ + fName = name ? name : "MHMcEnergyImpact"; + fTitle = title ? title : "Impact (radius) vs Energy distribution"; + + // - we initialize the histogram + // - we have to give diferent names for the diferent histograms because + // root don't allow us to have diferent histograms with the same name + + fHist.SetName(fName); + fHist.SetTitle(fTitle); + + fHist.SetDirectory(NULL); + + fHist.SetXTitle("E [GeV]"); + fHist.SetYTitle("r [m]"); + fHist.SetZTitle("N"); + + MBinning binsx; + binsx.SetEdgesLog(10, 1, 100000); // [GeV] + + MBinning binsy; + binsy.SetEdges(9, 0, 450); // [m] + + SetBinning(&fHist, &binsx, &binsy); + +} + +void MHMcEnergyImpact::SetName(const char *name) +{ + fName = name; + fHist.SetName(name); + fHist.SetDirectory(NULL); +} + +void MHMcEnergyImpact::SetTitle(const char *title) +{ + fTitle = title; + fHist.SetName(title); +} + +//------------------------------------------------------------------------- +// +// Defualt Destructor +// +MHMcEnergyImpact::~MHMcEnergyImpact() +{ +} + + +Bool_t MHMcEnergyImpact::SetupFill(const MParList *pList) +{ + const MBinning *binsx = (MBinning*)pList->FindObject("BinningEnergy", "MBinning"); + const MBinning *binsy = (MBinning*)pList->FindObject("BinningImpact", "MBinning"); + + if (!binsx || !binsy) + return kTRUE; + + SetBinning(&fHist, binsx, binsy); + + return kTRUE; +} + +//-------------------------------------------------------------------------- +// +// Fill the histogram with the log10 of the energy for triggered events. +// +Bool_t MHMcEnergyImpact::Fill(const MParContainer *cont) +{ + const MMcEvt &mcevt = *(MMcEvt*)cont; + + const Float_t energy = mcevt.GetEnergy(); + const Float_t impact = mcevt.GetImpact()/100.; + + fHist.Fill(energy, impact); + + return kTRUE; +} + +// ------------------------------------------------------------------------ +// +// Drawing function. It creates its own canvas. +// +void MHMcEnergyImpact::Draw(Option_t *option) +{ + if (!gPad) + MH::MakeDefCanvas(&fHist); + + gPad->SetLogx(); + + fHist.Draw(option); + + gPad->Modified(); + gPad->Update(); +} + +TObject *MHMcEnergyImpact::DrawClone(Option_t *option) const +{ + TCanvas *c = MH::MakeDefCanvas(&fHist); + + c->SetLogx(); + + // + // This is necessary to get the expected bahviour of DrawClone + // + gROOT->SetSelectedPad(NULL); + + ((TH2D&)fHist).DrawCopy(option); + + c->Modified(); + c->Update(); + + return c; +} + Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyImpact.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyImpact.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyImpact.h (revision 1974) @@ -0,0 +1,39 @@ +#ifndef MARS_MHMcEnergyImpact +#define MARS_MHMcEnergyImpact + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef ROOT_TH2 +#include +#endif + +class MMcEvt; + +class MHMcEnergyImpact : public MH +{ +private: + TH2D fHist; // histogram with the logarith of the energy + +public: + MHMcEnergyImpact(const char *name=NULL, const char *title=NULL); + ~MHMcEnergyImpact(); + + void SetName(const char *name); + void SetTitle(const char *title); + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *pcont=NULL); + + const TH2D *GetHist() { return &fHist; } + const TH2D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t* option = ""); + TObject *DrawClone(Option_t* option = "") const; + + ClassDef(MHMcEnergyImpact, 1) // Histogram container for 2D histogram in Energy and Impact +}; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.cc (revision 1974) @@ -0,0 +1,226 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek 4/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHMcEnergyMigration // +// // +// calculates the migration matrix E-est vs. E-true // +// for different bins in Theta // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHMcEnergyMigration.h" + +#include + +#include "MMcEvt.hxx" + +#include "MEnergyEst.h" +#include "MBinning.h" +#include "MHillasSrc.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHMcEnergyMigration); + + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHMcEnergyMigration::MHMcEnergyMigration(const char *name, const char *title) + : fHist() +{ + // + // set the name and title of this object + // + fName = name ? name : "MHMcEnergyMigration"; + fTitle = title ? title : "3-D histogram E-true E-est Theta"; + + fHist.SetDirectory(NULL); + + fHist.SetTitle("3D-plot E-true E-est Theta"); + fHist.SetXTitle("E_{true} [GeV]"); + fHist.SetYTitle("E_{est} [GeV]"); + fHist.SetZTitle("\\Theta [\\circ]"); +} + +// -------------------------------------------------------------------------- +// +// Set the binnings and prepare the filling of the histograms +// +Bool_t MHMcEnergyMigration::SetupFill(const MParList *plist) +{ + fEnergy = (MEnergyEst*)plist->FindObject("MEnergyEst"); + if (!fEnergy) + { + *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + const MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); + const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); + if (!binsenergy || !binstheta) + { + *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; + return kFALSE; + } + + SetBinning(&fHist, binsenergy, binsenergy, binstheta); + + fHist.Sumw2(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Draw a copy of the histogram +// +TObject *MHMcEnergyMigration::DrawClone(Option_t *opt) const +{ + TCanvas &c = *MakeDefCanvas("EnergyMigration", "E-est vs. E-true"); + + c.Divide(2, 2); + + gROOT->SetSelectedPad(NULL); + + TH1 *h; + + c.cd(1); + h = ((TH3*)(&fHist))->Project3D("expro"); + + h->SetTitle("Distribution of E-true"); + h->SetXTitle("E_{true} [GeV]"); + h->SetYTitle("Counts"); + + h->Draw(opt); + h->SetBit(kCanDelete); + gPad->SetLogx(); + + + c.cd(2); + h = ((TH3*)(&fHist))->Project3D("eypro"); + + h->SetTitle("Distribution of E-est"); + h->SetXTitle("E_{est} [GeV]"); + h->SetYTitle("Counts"); + + h->Draw(opt); + h->SetBit(kCanDelete); + gPad->SetLogx(); + + c.cd(3); + h = ((TH3*)(&fHist))->Project3D("ezpro"); + + h->SetTitle("Distribution of Theta"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetYTitle("Counts"); + + h->Draw(opt); + h->SetBit(kCanDelete); + + + c.cd(4); + ((TH3*)(&fHist))->DrawCopy(opt); + + c.Modified(); + c.Update(); + + return &c; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHMcEnergyMigration::Draw(Option_t *opt) +{ + if (!gPad) + MakeDefCanvas("EnergyMigration", "E-est vs. E-true"); + + gPad->Divide(2,2); + + TH1 *h; + + gPad->cd(1); + h = fHist.Project3D("ex_pro"); + + h->SetTitle("Distribution of E-true"); + h->SetXTitle("E_{true} [GeV]"); + + h->Draw(opt); + h->SetBit(kCanDelete); + gPad->SetLogx(); + + + gPad->cd(2); + h = fHist.Project3D("ey_pro"); + + h->SetTitle("Distribution of E-est"); + h->SetXTitle("E_{est} [GeV]"); + h->Draw(opt); + h->SetBit(kCanDelete); + gPad->SetLogx(); + + gPad->cd(3); + h = fHist.Project3D("ez_pro"); + + h->SetTitle("Distribution of Theta"); + h->SetXTitle("\\Theta [\\circ]"); + h->Draw(opt); + h->SetBit(kCanDelete); + + + gPad->cd(4); + fHist.DrawCopy(opt); + + gPad->Modified(); + gPad->Update(); +} + +// -------------------------------------------------------------------------- +// +// Fill the histogram +// +Bool_t MHMcEnergyMigration::Fill(const MParContainer *par) +{ + // MHillasSrc &hil = *(MHillasSrc*)par; + // fHist.Fill(fMcEvt->GetEnergy(), hil.GetSize()); + + // get E-true from fMcEvt and E-est from fEnergy + fHist.Fill(fMcEvt->GetEnergy(), fEnergy->GetEnergy(), fMcEvt->GetTheta()*kRad2Deg); + + return kTRUE; +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.h (revision 1974) @@ -0,0 +1,49 @@ +#ifndef MARS_MHMcEnergyMigration +#define MARS_MHMcEnergyMigration + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH3 +#include "TH3.h" +#endif + +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +class MMcEvt; +class MEnergyEst; +class MParList; + +class MHMcEnergyMigration : public MH +{ +private: + MMcEvt *fMcEvt; + MEnergyEst *fEnergy; + + TH3D fHist; + +public: + MHMcEnergyMigration(const char *name=NULL, const char *title=NULL); + + virtual Bool_t SetupFill(const MParList *pList); + virtual Bool_t Fill(const MParContainer *par); + + const TH3D *GetHist() { return &fHist; } + const TH3D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option="") const; + + ClassDef(MHMcEnergyMigration, 0) //3D-histogram E-true E-est Theta + +}; + +#endif + + + Index: /trunk/MagicSoft/Mars/mhistmc/MHMcIntRate.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcIntRate.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcIntRate.cc (revision 1974) @@ -0,0 +1,130 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Lopez 05/2001 +! Author(s): Thomas Bretz 05/2001 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMcIntRate +// +// This class holds the information (histogram and fit function) +// about the energy threshold for a particular trigger condition. +// +//////////////////////////////////////////////////////////////////////////// +#include "MHMcIntRate.h" + +#include + +#include + +#include "MH.h" +#include "MBinning.h" + +#include "MHMcDifRate.h" + +ClassImp(MHMcIntRate); + +// ------------------------------------------------------------------------- +// +// Default Constructor. +// +MHMcIntRate::MHMcIntRate(const char *name, const char *title) + : fHist() +{ + fName = name ? name : "MHMcIntRate"; + fTitle = title ? title : "Integral Trigger Rate"; + + // - we initialize the histogram + // - we have to give diferent names for the diferent histograms because + // root don't allow us to have diferent histograms with the same name + + fHist.SetName(fName); + fHist.SetTitle(fTitle); + + fHist.SetDirectory(NULL); + + fHist.SetXTitle("E [GeV]"); + fHist.SetYTitle("Rate [Hz]"); +} + +void MHMcIntRate::SetName(const char *name) +{ + fName = name; + fHist.SetName(name); + fHist.SetDirectory(NULL); +} + +void MHMcIntRate::SetTitle(const char *title) +{ + fTitle = title; + fHist.SetTitle(title); +} + +// ------------------------------------------------------------------------ +// +// Drawing function. It creates its own canvas. +// +void MHMcIntRate::Draw(Option_t *option) +{ + TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(""); + + pad->SetLogx(); + + fHist.Draw(option); + + pad->Modified(); + pad->Update(); +} + +// -------------------------------------------------------------------------- +// +// Calculate the IntRate and set the 'ReadyToSave' flag. +// The IntRate is calculated as the number of selected showers +// (eg. triggered ones) divided by the number of all showers. +// For error calculation Binomial errors are computed. +// +void MHMcIntRate::Calc(const MHMcDifRate &rate) +{ + /*const*/ TH1D &hist = (TH1D&)*rate.GetHist(); + const TAxis &axis = *hist.GetXaxis(); + + // + // Set the binning from the two axis of one of the two histograms + // + MH::SetBinning(&fHist, &axis); + + const Int_t nbinsx = axis.GetNbins(); + + for (Int_t i=1; i<=nbinsx; i++) + { + fHist.SetBinContent(i, hist.Integral(i, nbinsx, "width")); + fHist.SetBinError(i, hist.GetBinError(i)*axis.GetBinWidth(i)); + } + + fHist.SetEntries(hist.GetEntries()); + + SetReadyToSave(); +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcIntRate.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcIntRate.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcIntRate.h (revision 1974) @@ -0,0 +1,39 @@ +#ifndef MARS_MHMcIntRate +#define MARS_MHMcIntRate + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +#ifndef ROOT_TH1 +#include +#endif + +class MHMcDifRate; + +class MHMcIntRate : public MParContainer +{ +private: + + TH1D fHist; // histogram with the logarith of the energy + + // const MMcEvt *fMcEvt; //! container to fill histogram from + +public: + MHMcIntRate(const char *name=NULL, const char *title=NULL); + + void SetName(const char *name); + void SetTitle(const char *title); + + const TH1D *GetHist() { return &fHist; } + const TH1D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t* option = ""); + + void Calc(const MHMcDifRate &rate); + + ClassDef(MHMcIntRate, 1) // Histogram container for integral event rate +}; + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/MHMcRate.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcRate.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcRate.cc (revision 1974) @@ -0,0 +1,257 @@ +/* ======================================================================== *\ +! +! * +! * 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/2003 +! +! Explanations on the rate calculation can be found in +! chapter 7 of the following diploma thesis: +! http://www.pd.infn.it/magic/tesi2.ps.gz (in Italian) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +#include "MHMcRate.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHMcRate); + +void MHMcRate::Init(const char *name, const char *title) +{ + fName = name ? name : "MMcTriggerRate"; + fTitle = title ? title : "Task to calc the collection area "; + + fPartId=0; // Type of particle + + fEnergyMax=0.0; // Maximum Energy (TeV) + fEnergyMin=1000000.0; // Minimum Energy (TeV) + + fSolidAngle = -1.; // Solid angle within which incident directions + // are distributed + fThetaMax=0.0; // Maximum theta angle of run + fThetaMin=370.0; // Minimum theta angle of run + fPhiMax=0.0; // Maximum phi angle of run + fPhiMin=370.0; // Minimum phi angle of run + + fImpactMax=0.0; // Maximum impact parameter + fImpactMin=100000.0; // Minimum impact parameter + + fBackTrig=-1.0; // Number of triggers from background + fBackSim=-1.0; // Number of simulated showers for the background + + fTriggerRate= -1.0; // Trigger rate in Hz + fTriggerRateError= -1.0; // Estimated error for the trigger rate in Hz + + fMultiplicity = -1; // Multiplicity of the trigger condition + fMeanThreshold = -1.; // Discr. threshold of the trigger condition + fTriggerCondNum = 0; // Trigger condition number within camera file + +} + +// -------------------------------------------------------------------------- +// +// default constructor +// fills all member data with initial values +// +MHMcRate::MHMcRate(const char *name, const char *title) +{ + Init(name, title); + + fSpecIndex=0.0; // dn/dE = k * e^{- fSpecIndex} + fFlux0=-1.0; // dn/dE = fFlux0 * E^{-a} + + fShowerRate= -1.0; // Showers rate in Hz + fShowerRateError=0.0; // Estimated error of shower rate in Hz +} + +// -------------------------------------------------------------------------- +// +// overloaded constructor I +// fills all member data with initial values and sets the rate of +// incident showers to ShowRate +// +MHMcRate::MHMcRate(Float_t showrate, + const char *name, const char *title) +{ + Init(name, title); + + fSpecIndex=0.0; // dn/dE = k * e^{- fSpecIndex} + fFlux0=-1.0; // dn/dE = fFlux0 * E^{-a} + + fShowerRate= showrate; // Showers rate in Hz + fShowerRateError=sqrt(showrate); // Estimated error of shower rate in Hz +} + +// -------------------------------------------------------------------------- +// +// overloaded constructor I +// fills all member data with initial values and sets the +// spectral index and the initial flux to SpecIndex and Flux0 +// +MHMcRate::MHMcRate(Float_t specindex, Float_t flux0, + const char *name, const char *title) +{ + Init(name, title); + + fSpecIndex=specindex; // dn/dE = k * e^{- fSpecIndex} + fFlux0=flux0; // dn/dE = fFlux0 * E^{-a} + + fShowerRate= -1.0; + fShowerRateError=0.0; +} + +// -------------------------------------------------------------------------- +// +// set the particle that produces the showers in the athmosphere +// +void MHMcRate:: SetParticle(UShort_t part) +{ + fPartId=part; +} + +// -------------------------------------------------------------------------- +// +// Set the information about trigger due only to the Night Sky Background: +// +void MHMcRate::SetBackground (Float_t showers, Float_t triggers) +{ + fBackTrig=showers; // Number of triggers from background + fBackSim=triggers; // Number of simulated showers for the background +} + +// -------------------------------------------------------------------------- +// +// set the parameters to compute the incident rate +// +void MHMcRate:: SetFlux(Float_t flux0, Float_t specindx) +{ + fFlux0=flux0; + fSpecIndex=specindx; + +} + +// -------------------------------------------------------------------------- +// +// set the incident rate +// +void MHMcRate:: SetIncidentRate(Float_t showerrate) +{ + fShowerRate=showerrate; +} + +// -------------------------------------------------------------------------- +// +// update the limits for energy, theta, phi and impact parameter +// +void MHMcRate::UpdateBoundaries(Float_t energy, Float_t theta, + Float_t phi, Float_t impact) +{ + // It updates the limit values + + if (fThetaMaxtheta) fThetaMin=theta; + + if (fPhiMaxphi) fPhiMin=phi; + + if (fImpactMaximpact) fImpactMin=impact; + + if (fEnergyMaxenergy) fEnergyMin=energy; +} + +// -------------------------------------------------------------------------- +// +// compute the trigger rate and set the ReadyToSave bit +// +void MHMcRate::CalcRate(Float_t trig, Float_t anal, Float_t simu) +{ + // It computes the trigger rate + + // First one computes the rate of incident showers. + const Double_t specidx = 1.0-fSpecIndex; + + const Double_t epowmax = pow(fEnergyMax, specidx); + const Double_t epowmin = pow(fEnergyMin, specidx); + + if (fShowerRate <= 0) + fShowerRate = fFlux0/specidx*(epowmax-epowmin); + + if (fSolidAngle < 0.) + fSolidAngle = (fPhiMax-fPhiMin)*(cos(fThetaMin)-cos(fThetaMax)); + + if (fPartId!=1) + fShowerRate *= fSolidAngle; + + fShowerRate *= TMath::Pi()*(fImpactMax/100.0*fImpactMax/100.0 - + fImpactMin/100.0*fImpactMin/100.0); + + fShowerRateError = sqrt(fShowerRate); + + // The simulated trigger time in the camera program is 160 ns: + // 9/10/2002, AM: Fixed error below in calculation of "anal2" + // ( added factor fShowerRate/simu ) + + const Double_t anal2 = 1.0-fShowerRate*(anal/simu)*160.0e-9; + const Double_t back2 = fBackSim*160.0e-9; + + // Then the trigger rate and its error is evaluated + if(fBackTrig<0){ + fTriggerRateError = sqrt((trig*fShowerRate*fShowerRate/(simu*simu)) + + (anal2*anal2*1/(fBackSim*back2*back2))); + fBackTrig=0; + } + else + fTriggerRateError = sqrt((trig*fShowerRate*fShowerRate/(simu*simu)) + + (anal2*anal2*fBackTrig/(back2*back2))); + + fTriggerRate = trig*fShowerRate/simu + anal2*fBackTrig/back2; + + SetReadyToSave(); +} + +// -------------------------------------------------------------------------- +// +// print the trigger rate +// +void MHMcRate::Print(Option_t *) const +{ + *fLog << all << "Incident rate " << fShowerRate << " Hz " << endl; + *fLog << "Multiplicity: " << fMultiplicity << ", Discr. threshold: " << fMeanThreshold << endl; + *fLog << "Trigger Rate " << fTriggerRate << " +- " << fTriggerRateError << " Hz" << endl; +} + +// -------------------------------------------------------------------------- +// +// draw the trigger rate +// +void MHMcRate::Draw(Option_t *) +{ + *fLog << all << dbginf << " - MHMcRate::Draw: To be iplemented" << endl; +} + +TObject *MHMcRate::DrawClone(Option_t *) const +{ + *fLog << all << dbginf << " - MHMcRate::DrawClone: To be iplemented" << endl; + return NULL; +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcRate.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcRate.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcRate.h (revision 1974) @@ -0,0 +1,103 @@ +#ifndef MARS_MHMcRate +#define MARS_MHMcRate + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MHMcRate : public MParContainer +{ + +private: + UShort_t fPartId; // Type of particle + + Float_t fEnergyMax; // Maximum Energy in TeV + Float_t fEnergyMin; // Minimum Energy in TeV + + Float_t fThetaMax; // Maximum theta angle of run + Float_t fThetaMin; // Minimum theta angle of run + Float_t fPhiMax; // Maximum phi angle of run + Float_t fPhiMin; // Minimum phi angle of run + + Float_t fSolidAngle; // Solid angle within which incident directions + // are distributed (sr) + + Float_t fImpactMax; // Maximum impact parameter (cm) + Float_t fImpactMin; // Minimum impact parameter (cm) + + Float_t fBackTrig; // Number of triggers from background + Float_t fBackSim; // Number of simulated showers for the background + + Float_t fSpecIndex; // dn/dE = k * e^{- fSpecIndex} + Float_t fFlux0; // dn/dE = fFlux0 * E^{-a} + + Float_t fShowerRate; // Showers rate in Hz + Float_t fShowerRateError; // Estimated error of shower rate in Hz + + Float_t fTriggerRate; // Trigger rate in Hz + Float_t fTriggerRateError; // Estimated error for the trigger rate in Hz + + Float_t fMeanThreshold; // Mean discriminator threshold (mV) of trigger + // pixels. + + Short_t fMultiplicity; // L1 trigger multiplicity. + + Short_t fTriggerCondNum; // Trigger condition number, for the case of + // running over camra files containing several. + + void Init(const char *name, const char *title); + +public: + + MHMcRate(const char *name=NULL, const char *title=NULL); + MHMcRate(Float_t showrate, + const char *name=NULL, const char *title=NULL); + MHMcRate(Float_t specindex, Float_t flux0, + const char *name=NULL, const char *title=NULL); + + void SetParticle(UShort_t part); + void SetBackground(Float_t showers, Float_t triggers); + void SetFlux(Float_t flux0, Float_t specindx); + void SetIncidentRate(Float_t showerrate); + + void SetImpactMax(Float_t Impact) {fImpactMax=Impact;} + void SetImpactMin(Float_t Impact) {fImpactMin=Impact;} + + void SetThetaMax(Float_t Theta) {fThetaMax=Theta;} + void SetThetaMin(Float_t Theta) {fThetaMin=Theta;} + void SetPhiMax(Float_t Phi) {fPhiMax=Phi;} + void SetPhiMin(Float_t Phi) {fPhiMin=Phi;} + + void SetSolidAngle(Float_t Solid) {fSolidAngle=Solid;} + void SetEnergyMax(Float_t Energy) {fEnergyMax=Energy;} + void SetEnergyMin(Float_t Energy) {fEnergyMin=Energy;} + + void SetMultiplicity(Short_t nMul) {fMultiplicity = nMul;} + void SetMeanThreshold(Float_t thresh) {fMeanThreshold = thresh;} + + void SetTriggerCondNum(Short_t num) {fTriggerCondNum = num;} + + void UpdateBoundaries(Float_t energy, Float_t theta, Float_t phi, Float_t impact); + + Float_t GetTriggerRate() {return fTriggerRate;} + Float_t GetTriggerRateError() {return fTriggerRateError;} + + Short_t GetMultiplicity() {return fMultiplicity;} + Float_t GetMeanThreshold() {return fMeanThreshold;} + Short_t GetTriggerCondNum() {return fTriggerCondNum;} + + + void CalcRate(Float_t trig, Float_t anal, Float_t simu); + + void Print(Option_t *o=NULL) const; + + void Draw(Option_t *o=NULL); + TObject *DrawClone(Option_t *o=NULL) const; + + ClassDef(MHMcRate, 1) // Data Container to calculate trigger rate +}; + +#endif + + + Index: /trunk/MagicSoft/Mars/mhistmc/MHMcTriggerLvl2.cc =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcTriggerLvl2.cc (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcTriggerLvl2.cc (revision 1974) @@ -0,0 +1,311 @@ +/* ======================================================================== *\ +! +! * +! * 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): Nicola Galante, 2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMcTriggerLvl2 +// +// This class contains histograms for the Trigger Level2 image parameters +// +///////////////////////////////////////////////////////////////////////////// + +#include "MHMcTriggerLvl2.h" + +#include + +#include +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MMcTriggerLvl2.h" +#include "MGeomCam.h" +#include "MBinning.h" + +/* + Please, DON'T USE IFDEFS IN SUCH A CONTEXT, Thomas. + -------------------------------------------------- +#ifndef COLOR_LINELPS +#define COLOR_LINELPS Int_t colorlps = 1 +COLOR_LINELPS; +#endif + +#ifndef COLOR_LINESBC +#define COLOR_LINESBC Int_t colorsbc = 1 +COLOR_LINESBC; +#endif + +#ifndef COLOR_LINEPS +#define COLOR_LINEPS Int_t colorps = 1 +COLOR_LINEPS; +#endif +*/ + +/* Use this insteadif you want to have some value which is the same for all + your instances of MHMcTriggerLvl2, if not define the values in the constructor + and remove the 'static' + */ + +Int_t MHMcTriggerLvl2::fColorLps = 1; +Int_t MHMcTriggerLvl2::fColorSbc = 1; +Int_t MHMcTriggerLvl2::fColorPs = 1; + +ClassImp(MHMcTriggerLvl2); + +// -------------------------------------------------------------------------- +// +// Setup three histograms for fLutPseudoSize, fPseudoSize, fSizeBiggerCell +// +MHMcTriggerLvl2::MHMcTriggerLvl2(const char *name, const char *title) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHMcTriggerLvl2"; + fTitle = title ? title : "Trigger L2 image parameters"; + + hfLutPseudoSize = new TH1F("hfLutPseudoSize", "number of compact pixels in one lut", 13, 0, 12); + hfPseudoSize = new TH1F("hfPseudoSize", "Multiplicity of the cluster identified by the L2T", 13, 0, 12); + hfSizeBiggerCell = new TH1F("hfSizeBiggerCell", "Number of fired pixel in bigger cell", 37, 0, 36); + + hfLutPseudoSizeNorm = new TH1F("hfLutPseudoSizeNorm", "Normalized Number of compact pixels in one lut", 13, 0, 12); + hfPseudoSizeNorm = new TH1F("hfPseudoSizeNorm", "Normalized Multiplicity of the cluster identified by the L2T", 13, 0, 12); + hfSizeBiggerCellNorm = new TH1F("hfSizeBiggerCellNorm", "Normalized Number of fired pixel in bigger cell", 37, 0, 36); + + hfLutPseudoSize->SetDirectory(NULL); + hfLutPseudoSizeNorm->SetDirectory(NULL); + hfPseudoSize->SetDirectory(NULL); + hfPseudoSizeNorm->SetDirectory(NULL); + hfSizeBiggerCell->SetDirectory(NULL); + hfSizeBiggerCellNorm->SetDirectory(NULL); + + hfLutPseudoSize->SetXTitle("Number of Pixels"); + hfLutPseudoSizeNorm->SetXTitle("Number of Pixels"); + hfPseudoSize->SetXTitle("Number of Pixels"); + hfPseudoSizeNorm->SetXTitle("Number of Pixels"); + hfSizeBiggerCell->SetXTitle("Number of Pixels"); + hfSizeBiggerCellNorm->SetXTitle("Number of Pixels"); + + hfLutPseudoSize->SetYTitle("Counts"); + hfLutPseudoSizeNorm->SetYTitle("Counts/Total Counts"); + hfPseudoSize->SetYTitle("Counts"); + hfPseudoSizeNorm->SetYTitle("Counts/Total Counts"); + hfSizeBiggerCell->SetYTitle("Counts"); + hfSizeBiggerCellNorm->SetYTitle("Counts/Total Counts"); + + f1 = new TF1("FNorm", "1", -1, 40); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms +// +MHMcTriggerLvl2::~MHMcTriggerLvl2() +{ + delete hfLutPseudoSize; + delete hfPseudoSize; + delete hfSizeBiggerCell; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms with data from a MMcTriggerLvl2-Container. +// Be careful: Only call this with an object of type MMcTriggerLvl2 +// +Bool_t MHMcTriggerLvl2::Fill(const MParContainer *par) +{ + const MMcTriggerLvl2 &h = *(MMcTriggerLvl2 *)par; + + hfLutPseudoSize->Fill(h.GetLutPseudoSize()); + hfPseudoSize->Fill(h.GetPseudoSize()); + hfSizeBiggerCell->Fill(h.GetSizeBiggerCell()); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// This is the private function member which draw a clone of a histogram. +// This method is called by the DrawClone method. +// +TObject *MHMcTriggerLvl2::DrawHist(TH1 &hist, TH1 &histNorm, const TString &canvasname, Int_t &col) const +{ + col++; + + TCanvas *c = (TCanvas*)gROOT->FindObject(canvasname); + + Bool_t same = kTRUE; + if (!c) + { + c = MakeDefCanvas(canvasname,canvasname, 800, 610); + c->Divide(2,1); + same = kFALSE; + } + + c->cd(1); + hist.SetLineColor(col); + hist.DrawCopy(same?"same":""); + + c->cd(2); + histNorm.SetLineColor(col); + histNorm.DrawCopy(same?"same":""); + + return c; +} + + +// -------------------------------------------------------------------------- +// +// Draw a clone of a data member histogram. So that the object can be deleted +// and the histogram is still visible in the canvas. +// The cloned object are deleted together with the canvas if the canvas is +// destroyed. If you want to handle dostroying the canvas you can get a +// pointer to it from this function +// Possible options are: +// "lps" for the fLutPseudoSize histogram; +// "sbc" for the fSizeBiggerCell histogram; +// "ps" for the fPseudoSize histogram; +// default: "ps" +// +TObject *MHMcTriggerLvl2::DrawClone(Option_t *opt) const +{ + TString str(opt); + + if (str.IsNull()) + str = "ps"; + + if (!str.Contains("lps", TString::kIgnoreCase) && + !str.Contains("sbc", TString::kIgnoreCase) && + !str.Contains("ps", TString::kIgnoreCase)) + { + *fLog << "ARGH!@! Possible options are \"lps\", \"sbc\", \"ps\" or NULL!" <Draw(); + pad2->Draw(); + pad3->Draw(); + + pad1->cd(); + hfLutPseudoSize->Draw(); + + pad2->cd(); + hfSizeBiggerCell->Draw(); + + pad3->cd(); + hfPseudoSize->Draw(); +} + + + +// -------------------------------------------------------------------------- +// +// Return the histogram by its name +// +TH1 *MHMcTriggerLvl2::GetHistByName(const TString name) +{ + if (name.Contains("hfLutPseudoSize", TString::kIgnoreCase)) + return hfLutPseudoSize; + if (name.Contains("hfSizeBiggerCell", TString::kIgnoreCase)) + return hfSizeBiggerCell; + if (name.Contains("hfPseudoSize", TString::kIgnoreCase)) + return hfPseudoSize; + + if (name.Contains("hfLutPseudoSizeNorm", TString::kIgnoreCase)) + return hfLutPseudoSizeNorm; + if (name.Contains("hfSizeBiggerCellNorm", TString::kIgnoreCase)) + return hfSizeBiggerCellNorm; + if (name.Contains("hfPseudoSizeNorm", TString::kIgnoreCase)) + return hfPseudoSizeNorm; + + return NULL; +} + + + +// -------------------------------------------------------------------------- +// +// Normalize the histogram on its integral, i.e. normalize the +// values of the histogram on the statistics. This method creates +// a copy (histNorm) of the histogram to normalize (hist) and normalize +// the copy (histNorm) on its integral. It returns histNorm. +// +TH1 *MHMcTriggerLvl2::NormalizeHist(TH1 *histNorm, TH1 *hist) const +{ + if (histNorm == hist){ + *fLog << "ARGH!@! You cannot pass the same histogram into each argument!" << endl; + return NULL; + } + + if ((hist == NULL) || (hist->Integral() == (Stat_t)0)){ + *fLog << "ARGH!@! You are trying to normalize the histogram to 0!" << endl; + return NULL; + } + + histNorm->Reset("ICE"); + histNorm->Add(hist, 1); + histNorm->Divide(f1, (Double_t)(hist->Integral())); + + return histNorm; +} Index: /trunk/MagicSoft/Mars/mhistmc/MHMcTriggerLvl2.h =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/MHMcTriggerLvl2.h (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/MHMcTriggerLvl2.h (revision 1974) @@ -0,0 +1,61 @@ +#ifndef MARS_MHMcTriggerLvl2 +#define MARS_MHMcTriggerLvl2 + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1F; +class TF1; +class MMcTriggerLvl2; + +class MHMcTriggerLvl2 : public MH +{ + public: + + +private: + + TH1F *hfLutPseudoSize; // Histogram of fLutPseudoSize + TH1F *hfLutPseudoSizeNorm; // Histogram of fLutPseudoSize normalized on integral of distribution + TH1F *hfPseudoSize; // Histogram of fPseudoSize + TH1F *hfPseudoSizeNorm; // Histogram of fPseudoSize normalized on integral of distribution + TH1F *hfSizeBiggerCell; // Histogram of fSizeBiggerCell + TH1F *hfSizeBiggerCellNorm; // Histogram of fSizeBiggerCell normalized on integral of distribution + + TF1* f1; // Function used to normalize histograms + + static Int_t fColorLps; + static Int_t fColorSbc; + static Int_t fColorPs; + +public: + MHMcTriggerLvl2(const char *name=NULL, const char *title=NULL); + ~MHMcTriggerLvl2(); + + Bool_t Fill(const MParContainer *par); + + TH1 *GetHistByName(const TString name); + + TH1F *GetHistfLutPseudoSize() const { return hfLutPseudoSize; } + TH1F *GetHistfLutPseudoSizeNorm() const { return hfLutPseudoSizeNorm; } + TH1F *GetHistfPseudoSize() const { return hfPseudoSize; } + TH1F *GetHistfPseudoSizeNorm() const { return hfPseudoSizeNorm; } + TH1F *GetHistfSizeBiggerCell() const { return hfSizeBiggerCell; } + TH1F *GetHistfSizeBiggerCellNorm() const { return hfSizeBiggerCellNorm; } + + void Draw(Option_t *opt=NULL); + TObject *DrawClone(Option_t *opt=NULL) const; + + TH1 *NormalizeHist(TH1 *histNorm, TH1 *hist) const; + + ClassDef(MHMcTriggerLvl2, 1) // Container which holds histograms for the Trigger Level2 image parameters + +private: + + TObject *DrawHist(TH1 &hist, TH1 &histNorm, const TString &canvasname, Int_t &colore) const; + +}; + + +#endif Index: /trunk/MagicSoft/Mars/mhistmc/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mhistmc/Makefile (revision 1974) +++ /trunk/MagicSoft/Mars/mhistmc/Makefile (revision 1974) @@ -0,0 +1,65 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +# +# Handling name of the Root Dictionary Files +# +CINT = HistMc + +# +# Library name to creatre +# +LIB = mhistmc.a + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mhist -I../mmc + +#------------------------------------------------------------------------------ + +.SUFFIXES: .c .cc .cxx .h .hxx .o + +SRCFILES = MHMcEnergy.cc \ + MHMcEfficiency.cc \ + MHMcEfficiencyImpact.cc \ + MHMcEfficiencyEnergy.cc \ + MHMcEnergyImpact.cc \ + MHMcRate.cc \ + MHMcCollectionArea.cc \ + MHMcIntRate.cc \ + MHMcDifRate.cc \ + MHMcEnergyMigration.cc \ + MHMcCT1CollectionArea.cc \ + MHMcTriggerLvl2.cc + +SRCS = $(SRCFILES) +HEADERS = $(SRCFILES:.cc=.h) +OBJS = $(SRCFILES:.cc=.o) + +############################################################ + +all: $(LIB) + +include ../Makefile.rules + +clean: rmcint rmobjs rmcore rmlib + +mrproper: clean rmbak + +# @endcode + + + + + + + +