Index: trunk/MagicSoft/Mars/mhist/HistLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 1414) +++ trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 1415) @@ -18,8 +18,9 @@ #pragma link C++ class MHEnergyTime+; #pragma link C++ class MHEnergyTheta+; -#pragma link C++ class MHAlphaEnergyTime; -#pragma link C++ class MHAlphaEnergyTheta; -#pragma link C++ class MHEffOnTimeTime+; -#pragma link C++ class MHEffOnTimeTheta+; +#pragma link C++ class MHAlphaEnergyTime+; +#pragma link C++ class MHAlphaEnergyTheta+; +#pragma link C++ class MHGamma+; +#pragma link C++ class MHEffOnTime+; + #pragma link C++ class MHTimeDiffTime+; #pragma link C++ class MHTimeDiffTheta+; @@ -40,7 +41,6 @@ #pragma link C++ class MHMcCollectionArea+; #pragma link C++ class MHMcEnergyMigration+; +#pragma link C++ class MHFlux; -#pragma link C++ class MHCompProb+; -#pragma link C++ class MHHadroness+; #endif Index: trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc (revision 1414) +++ trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc (revision 1415) @@ -66,5 +66,5 @@ fHist.SetTitle("3D-plot of alpha, E-est, Theta"); fHist.SetXTitle("\\alpha [\\circ]"); - fHist.SetYTitle("E_{est} [GeV]"); + fHist.SetYTitle("E-est [GeV] "); fHist.SetZTitle("\\Theta [\\circ]"); } @@ -79,5 +79,5 @@ if (!fEnergy) { - *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl; + *fLog << err << dbginf << "MHAlphaEnergyTheta : MEnergyEst not found... aborting." << endl; return kFALSE; } @@ -86,21 +86,23 @@ if (!fMcEvt) { - *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + *fLog << err << dbginf << "MHAlphaEnergyTheta : MMcEvt not found... aborting." << endl; return kFALSE; } MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); - MBinning* binsalpha = (MBinning*)plist->FindObject("BinningAlpha"); + MBinning* binsalphaflux = (MBinning*)plist->FindObject("BinningAlphaFlux"); MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); - if (!binsenergy || !binsalpha || !binstheta) + if (!binsenergy || !binsalphaflux || !binstheta) { - *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; + *fLog << err << dbginf << "MHAlphaEnergyTheta : At least one MBinning not found... aborting." << endl; return kFALSE; } - SetBinning(&fHist, binsalpha, binsenergy, binstheta); + SetBinning(&fHist, binsalphaflux, binsenergy, binstheta); fHist.Sumw2(); + fHist.Sumw2(); + return kTRUE; } @@ -114,5 +116,5 @@ MHillasSrc &hil = *(MHillasSrc*)par; - fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(), fMcEvt->GetTheta()*kRad2Deg); + fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(), fMcEvt->GetTheta()*kRad2Deg); return kTRUE; @@ -146,5 +148,5 @@ h->SetTitle("Distribution of E-est [GeV]"); - h->SetXTitle("E_{est} [GeV]"); + h->SetXTitle("E-est [GeV] "); h->SetYTitle("Counts"); @@ -223,134 +225,2 @@ } -// -------------------------------------------------------------------------- -// -// Calculate the histogram as the difference of two histograms : -// fHist(gamma) = h1(source) - h2(antisource) -// -void MHAlphaEnergyTheta::Subtract(const TH3D *h1, const TH3D *h2) -{ - // MH::SetBinning(&fHist, (TH1*)h1); - - // fHist.Sumw2(); - fHist.Add((TH1*)h1, (TH1*)h2, 1, -1); // Root: FIXME -} - - -// -------------------------------------------------------------------------- -// -// Integrate fHist(gamma) in the alpha range (lo, up) -// -TH2D *MHAlphaEnergyTheta::GetAlphaProjection(Axis_t lo, Axis_t up) -{ - if (up < lo) - { - *fLog << err << fName << ": Alpha projection not possible: lo=" << lo << " up=" << up << endl; - return NULL; - } - - TAxis &axe = *fHist.GetXaxis(); - - Int_t ilo = axe.FindFixBin(lo); - Int_t iup = axe.FindFixBin(up); - - const Double_t epslo1 = lo-axe.GetBinLowEdge(ilo); - const Double_t epslo2 = axe.GetBinUpEdge(ilo)-lo; - - const Double_t epsup1 = up-axe.GetBinLowEdge(iup); - const Double_t epsup2 = axe.GetBinUpEdge(iup)-up; - - const Double_t epslo = epslo1epslo2) - ilo++; - - if (epsup10.01*axe.GetBinWidth(ilo) || epsup>0.01*axe.GetBinWidth(iup)) - { - *fLog << err << fName << ": binning is not adequate for the requested projection:" << endl; - *fLog << "Please specify a lower or upper limit which is not more than 1% away from a bin edge" << endl; - *fLog << " epslo = " << epslo << endl; - *fLog << " epsup = " << epsup << endl; - *fLog << " dwl = " << axe.GetBinWidth(ilo) << endl; - *fLog << " dwu = " << axe.GetBinWidth(iup) << endl; - return NULL; - } - - axe.SetRange(ilo, iup); - - TH2D &h2D = *(TH2D *)fHist.Project3D("ezypro"); - - h2D.SetTitle("2D-plot of Theta vs. E-est"); - h2D.SetXTitle("E-est [GeV] "); - h2D.SetYTitle("\\Theta [\\circ]"); - - return &h2D; -} - -// -------------------------------------------------------------------------- -// -// Draw the integrated histogram -// -TH2D *MHAlphaEnergyTheta::DrawAlphaProjection(Axis_t lo, Axis_t up, Option_t *opt) -{ - TH2D *h2D = GetAlphaProjection(lo, up); - - if (!h2D) - return NULL; - - char txt[100]; - sprintf(txt, "No.of Gammas vs. E-est and Theta (%.1f < alpha < %.1f deg)", lo, up); - - // TCanvas *c = MakeDefCanvas("AlphaEnergyTheta", "2D histogram of gamma signal in energy and theta"); - TCanvas &c = *MakeDefCanvas("AlphaEnergyTheta", txt); - - c.Divide(2, 2); - - gROOT->SetSelectedPad(NULL); - - TH1 *h; - - c.cd(1); - h = h2D->ProjectionX("Eest", -1, 9999, "E"); - h->SetTitle("Distribution of E-est [GeV]"); - h->SetXTitle("E-est [GeV] "); - h->SetYTitle("Counts"); - - h->Draw(opt); - h->SetBit(kCanDelete); - gPad->SetLogx(); - - c.cd(2); - h = h2D->ProjectionY("theta", -1, 9999, "E"); - h->SetTitle("Distribution of \\Theta [\\circ]"); - h->SetXTitle("\\Theta [\\circ]"); - h->SetYTitle("Counts"); - - h->Draw(opt); - h->SetBit(kCanDelete); - - c.cd(3); - - h2D->DrawCopy(opt); - gPad->SetLogx(); - - c.Modified(); - c.Update(); - - return h2D; -} - - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h (revision 1414) +++ trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h (revision 1415) @@ -42,13 +42,4 @@ TObject *DrawClone(Option_t *option="") const; - void Subtract(const TH3D *h1, const TH3D *h2); - void Subtract(const MHAlphaEnergyTheta *h1, const MHAlphaEnergyTheta *h2) - { - Subtract(h1->GetHist(), h2->GetHist()); - } - - TH2D *DrawAlphaProjection(Axis_t lo, Axis_t up, Option_t *opt=""); - TH2D *GetAlphaProjection(Axis_t lo, Axis_t up); - ClassDef(MHAlphaEnergyTheta, 1) //3D-histogram in alpha, Energy and theta }; Index: trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.cc (revision 1414) +++ trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.cc (revision 1415) @@ -68,5 +68,5 @@ fHist.SetTitle("3D-plot of alpha, E-est, time"); fHist.SetXTitle("\\alpha [\\circ]"); - fHist.SetYTitle("E_{est} [GeV]"); + fHist.SetYTitle("E-est [GeV] "); fHist.SetZTitle("time [s]"); } @@ -81,5 +81,5 @@ if (!fEnergy) { - *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl; + *fLog << err << dbginf << "MHAlphaEnergyTime : MEnergyEst not found... aborting." << endl; return kFALSE; } @@ -88,18 +88,18 @@ if (!fTime) { - *fLog << err << dbginf << "MTime not found... aborting." << endl; + *fLog << err << dbginf << "MHAlphaEnergyTime : MTime not found... aborting." << endl; return kFALSE; } MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); - MBinning* binsalpha = (MBinning*)plist->FindObject("BinningAlpha"); + MBinning* binsalphaflux = (MBinning*)plist->FindObject("BinningAlphaFlux"); MBinning* binstime = (MBinning*)plist->FindObject("BinningTime"); - if (!binsenergy || !binsalpha || !binstime) + if (!binsenergy || !binsalphaflux || !binstime) { - *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; + *fLog << err << dbginf << "MHAlphaEnergyTime : At least one MBinning not found... aborting." << endl; return kFALSE; } - SetBinning(&fHist, binsalpha, binsenergy, binstime); + SetBinning(&fHist, binsalphaflux, binsenergy, binstime); fHist.Sumw2(); @@ -116,5 +116,5 @@ MHillasSrc &hil = *(MHillasSrc*)par; - fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(), 0.0001*fTime->GetTimeLo()); + fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(), 0.0001*fTime->GetTimeLo()); return kTRUE; } @@ -147,5 +147,5 @@ h->SetTitle("Distribution of E-est [GeV]"); - h->SetXTitle("E_{est} [GeV]"); + h->SetXTitle("E-est [GeV] "); h->SetYTitle("Counts"); @@ -226,122 +226,110 @@ } -// -------------------------------------------------------------------------- -// -// Calculate the histogram as the difference of two histograms : -// fHist(gamma) = h1(source) - h2(antisource) -// -void MHAlphaEnergyTime::Subtract(const TH3D *h1, const TH3D *h2) -{ - // MH::SetBinning(&fHist, (TH1*)h1); - - // fHist.Sumw2(); - fHist.Add((TH1*)h1, (TH1*)h2, 1, -1); // ROOT: FIXME! -} - -// -------------------------------------------------------------------------- -// -// Integrate fHist(gamma) in the alpha range (lo, up) -// -TH2D *MHAlphaEnergyTime::GetAlphaProjection(Axis_t lo, Axis_t up) -{ - if (up < lo) + +// -------------------------------------------------------------------------- +// +// Integrate fHist (Alpha,E-est,Time) over the Time to get +// fAlphaEest(Alpha,E-est) +// +TH2D *MHAlphaEnergyTime::IntegrateTime(const char *title, Bool_t draw) +{ + Int_t nzbins = fHist.GetNbinsZ(); + TAxis &axez = *fHist.GetZaxis(); + axez.SetRange(1,nzbins); + + TH2D &fAlphaEest = *(TH2D *)fHist.Project3D("exy"); + + fAlphaEest.SetTitle(title); + fAlphaEest.SetXTitle("E-est [GeV] "); + fAlphaEest.SetYTitle("\\alpha [ \\circ]"); + + if (draw == kTRUE) { - *fLog << err << fName << ": Alpha projection not possible: lo=" << lo << " up=" << up << endl; - return NULL; + TCanvas &c = *MakeDefCanvas(title, title); + + gROOT->SetSelectedPad(NULL); + + fAlphaEest.DrawCopy(); + gPad->SetLogx(); + + c.Modified(); + c.Update(); } - TAxis &axe = *fHist.GetXaxis(); - - Int_t ilo = axe.FindFixBin(lo); - Int_t iup = axe.FindFixBin(up); - - const Double_t epslo1 = lo-axe.GetBinLowEdge(ilo); - const Double_t epslo2 = axe.GetBinUpEdge(ilo)-lo; - - const Double_t epsup1 = up-axe.GetBinLowEdge(iup); - const Double_t epsup2 = axe.GetBinUpEdge(iup)-up; - - const Double_t epslo = epslo1epslo2) - ilo++; - - if (epsup10.01*axe.GetBinWidth(ilo) || epsup>0.01*axe.GetBinWidth(iup)) + return &fAlphaEest; +} + +// -------------------------------------------------------------------------- +// +// Integrate fHist (Alpha,E-est,Time) over E-est to get +// fAlphaTime(Alpha,Time) +// +TH2D *MHAlphaEnergyTime::IntegrateEest(const char *title, Bool_t draw) +{ + Int_t nybins = fHist.GetNbinsY(); + TAxis &axey = *fHist.GetYaxis(); + axey.SetRange(1,nybins); + + TH2D &fAlphaTime = *(TH2D *)fHist.Project3D("exz"); + + fAlphaTime.SetTitle(title); + fAlphaTime.SetXTitle("Time [s]"); + fAlphaTime.SetYTitle("\\alpha [ \\circ]"); + + if (draw == kTRUE) { - *fLog << err << fName << ": binning is not adequate for the requested projection:" << endl; - *fLog << "Please specify a lower or upper limit which is not more than 1% away from a bin edge" << endl; - *fLog << " epslo = " << epslo << endl; - *fLog << " epsup = " << epsup << endl; - *fLog << " dwl = " << axe.GetBinWidth(ilo) << endl; - *fLog << " dwu = " << axe.GetBinWidth(iup) << endl; - return NULL; + TCanvas &c = *MakeDefCanvas(title, title); + + gROOT->SetSelectedPad(NULL); + + fAlphaTime.DrawCopy(); + + c.Modified(); + c.Update(); } - axe.SetRange(ilo, iup); - - TH2D &h2D = *(TH2D *)fHist.Project3D("ezy"); - - h2D.SetTitle("2D-plot of time vs. E-est"); - h2D.SetXTitle("E-est [GeV] "); - h2D.SetYTitle("time [s]"); - - return &h2D; -} - -//--------------------------------------------------------- -// -// Draw the projected histogram -// -TH2D *MHAlphaEnergyTime::DrawAlphaProjection(Axis_t lo, Axis_t up, Option_t *opt) -{ - TH2D *h2D = GetAlphaProjection(lo, up); - - if (!h2D) - return NULL; - - char txt[100]; - sprintf(txt, "No.of Gammas vs. E-est and Time (%.1f < alpha < %.1f deg)", lo, up); - - // TCanvas *c = MakeDefCanvas("AlphaEnergyTime", "2D histogram of gamma signal in energy and time"); - TCanvas &c = *MakeDefCanvas("AlphaEnergyTime", txt); - - c.Divide(2, 2); - - gROOT->SetSelectedPad(NULL); - - TH1 *h; - - c.cd(1); - h = h2D->ProjectionX("Eest", -1, 9999, "E"); - h->SetTitle("Distribution of E-est [GeV]"); - h->SetXTitle("E-est [GeV] "); - h->SetYTitle("Counts"); - - h->Draw(opt); - h->SetBit(kCanDelete); - gPad->SetLogx(); - - c.cd(2); - h = h2D->ProjectionY("time", -1, 9999, "E"); - h->SetTitle("Distribution of time [s]"); - h->SetXTitle("time [s]"); - h->SetYTitle("Counts"); - - h->Draw(opt); - h->SetBit(kCanDelete); - - c.cd(3); - - h2D->DrawCopy(opt); - gPad->SetLogx(); - - c.Modified(); - c.Update(); - - return h2D; -} - + return &fAlphaTime; +} + +// -------------------------------------------------------------------------- +// +// Integrate fHist (Alpha,E-est,Time) over Eest and Time to get +// fAlpha(Alpha) +// +TH1D *MHAlphaEnergyTime::IntegrateEestTime(const char *title, Bool_t draw) +{ + Int_t nybins = fHist.GetNbinsY(); + TAxis &axey = *fHist.GetYaxis(); + axey.SetRange(1,nybins); + + Int_t nzbins = fHist.GetNbinsZ(); + TAxis &axez = *fHist.GetZaxis(); + axez.SetRange(1,nzbins); + + TH1D &fAlpha = *(TH1D *)fHist.Project3D("ex"); + + fAlpha.SetTitle(title); + fAlpha.SetXTitle("\\alpha [ \\circ]"); + fAlpha.SetYTitle("Counts"); + + if (draw == kTRUE) + { + TCanvas &c = *MakeDefCanvas(title, title); + + gROOT->SetSelectedPad(NULL); + + fAlpha.DrawCopy(); + + c.Modified(); + c.Update(); + } + + return &fAlpha; +} + + + + + + + Index: trunk/MagicSoft/Mars/mhist/MHFlux.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFlux.cc (revision 1415) +++ trunk/MagicSoft/Mars/mhist/MHFlux.cc (revision 1415) @@ -0,0 +1,705 @@ +/* ======================================================================== *\ +! +! * +! * 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 5/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHFlux // +// // +// calculates absolute photon fluxes // +// from the distributions of the estimated energy // +// for the different bins in some variable 'Var' // +// (Var = Theta or time) // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHFlux.h" + +#include + +#include +#include +#include + + +#include + +#include "MTime.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHFlux); + + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets the variable name (Theta or time) +// and the units for the variable +// +MHFlux::MHFlux(const TH2D &h2d, const Bool_t Draw, + const char *varname, const char *unit) + : fHOrig(), fHUnfold(), fHFlux() +{ + if (varname == NULL || unit == NULL) + { + *fLog << "MHFlux : varname or unit not defined" << endl; + } + + fVarname = varname; + fUnit = unit; + + TString strg(varname); + strg += unit; + + // char txt[100]; + + // original distribution of E-est for different bins + // of the variable (Theta or time) + // sprintf(txt, "gammas vs. E-est and %s",varname); + + TString strg1 = "no.of gammas vs. E-est and "; + strg1 += varname; + + ((TH2D&)h2d).Copy(fHOrig); + + fHOrig.SetName("E-est"); + fHOrig.SetTitle(strg1); + + fHOrig.SetDirectory(NULL); + fHOrig.SetXTitle("E-est [GeV] "); + fHOrig.SetYTitle(strg); + fHOrig.Sumw2(); + + SetBinning((TH2*)&fHOrig, (TH2*)&h2d); + + + // unfolded distribution of E-unfold for different bins + // of the variable (Theta or time) + // sprintf(txt, "gammas vs. E-unfold and %s",varname); + TString strg2 = "no.of gammas vs. E-unfold and "; + strg2 += varname; + + fHUnfold.SetName("E-unfolded"); + fHUnfold.SetTitle(strg2); + + fHUnfold.SetDirectory(NULL); + fHUnfold.SetXTitle("E-unfold [GeV] "); + fHUnfold.SetYTitle(strg); + fHUnfold.Sumw2(); + + SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig); + + + // absolute photon flux vs. E-unfold + // for different bins of the variable (Theta or time) + // + // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname); + TString strg3 = "gamma flux [1/(s m2 GeV) vs. E-unfold and "; + strg3 += varname; + + fHFlux.SetName("photon flux"); + fHFlux.SetTitle(strg3); + + fHFlux.SetDirectory(NULL); + fHFlux.SetXTitle("E-unfold [GeV] "); + fHFlux.SetYTitle(strg); + fHFlux.Sumw2(); + + SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold); + + + // copy fHOrig into fHUnfold in case no unfolding is done + const Int_t nEunf = fHUnfold.GetNbinsX(); + const Int_t nVar = fHUnfold.GetNbinsY(); + for (int m=1; m<=nEunf; m++) + { + for (int n=1; n<=nVar; n++) + { + Double_t cont = fHOrig.GetBinContent(m,n); + Double_t dcont = fHOrig.GetBinError(m,n); + fHUnfold.SetBinContent(m,n,cont); + fHUnfold.SetBinError(m,n,dcont); + } + } + //.............................................. + // draw the No.of photons vs. E-est + // for the individual bins of the variable Var + + if (Draw == kTRUE) + { + const Int_t nVar = fHOrig.GetNbinsY(); + + for (int n=1; n<=nVar; n++) + { + TString strg0("Orig-"); + strg0 += fVarname; + TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E"); + + char txt0[100]; + strg0 = fVarname; + strg0 += "-bin %d"; + sprintf(txt0, strg0, n); + + TString strg1("No.of photons vs. E-est for "); + strg1 += txt0; + new TCanvas(txt0,strg1); + // TCanvas &c = *MakeDefCanvas(txt0, strg1); + // gROOT->SetSelectedPad(NULL); + + gPad->SetLogx(); + + h.SetName(txt0); + h.SetTitle(strg1); + h.SetXTitle("E-est [GeV] "); + h.SetYTitle("No.of photons"); + h.DrawCopy(); + + // c.Modified(); + // c.Update(); + } + } + //........................ +} + +// ------------------------------------------------------------------------- +// +// Dummy Fill (has to be included because in base class MH Fill is set to 0); +// without the dummy Fill one gets the error message : +// +// Error: Can't call MHFlux::MHFlux(evttime,"time","[s]") in current scope +// FILE:macros/flux.C LINE:465 +// Possible candidates are... +// filename line:size busy function type and name (in MHFlux) +// filename line:size busy function type and name (in MH) +// filename line:size busy function type and name (in MParContainer) +// filename line:size busy function type and name (in TObject) +// +Bool_t MHFlux::Fill(const MParContainer *par) +{ + return kTRUE; +} + + +// ------------------------------------------------------------------------- +// +// Unfold the distribution in E-est +// +void MHFlux::Unfold(const Bool_t Draw) +{ + //.............................................. + // draw the No.of photons vs. E-unfold + // for the individual bins of the variable Var + + if (Draw == kTRUE) + { + const Int_t nVar = fHUnfold.GetNbinsY(); + + for (int n=1; n<=nVar; n++) + { + TString strg0("Unfold-"); + strg0 += fVarname; + TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E"); + + char txt0[100]; + strg0 = fVarname; + strg0 += "-bin %d"; + sprintf(txt0, strg0, n); + + TString strg1("No.of photons vs. E-unfold for "); + strg1 += txt0; + new TCanvas(txt0,strg1); + + // TCanvas &c = *MakeDefCanvas(txt0, strg1); + // gROOT->SetSelectedPad(NULL); + + gPad->SetLogx(); + + h.SetName(txt0); + h.SetTitle(strg1); + h.SetXTitle("E-unfold [GeV] "); + h.SetYTitle("No.of photons"); + h.DrawCopy(); + + // c.Modified(); + // c.Update(); + } + } + //........................ +} + + +// ------------------------------------------------------------------------- +// +// Calculate photon flux by dividing the distribution in Eunf (fHUnfold) by +// the width of the energy interval (deltaE) +// the effective ontime (*teff) +// and the effective collection area (*aeff) +// +void MHFlux::CalcFlux(const TH1D *teff, const TProfile *thetabar, + const TH2D *aeff, const Bool_t Draw) +{ + // Note that fHUnfold has bins in Eunf and Var + // *teff has bins in Var (the same bins in Var as fHUnfold) + // *thetabar has bins in Var (the same bins in Var as fHUnfold) + // *aeff has bins in Etru and Theta + // (where in general the binning in Etru is different + // from the binning in Eunf) + // The variable Var may be 'time' or 'Theta' + + // Draw = kTRUE means the differential flux vs E-unf should be drawn + // for the individual bins of the variable Var + + //.................................... + // define dummy histogram *aeff + ((TH1*)aeff)->Sumw2(); + MBinning binsetru("BinningEtru"); + binsetru.SetEdgesLog(10, 10, 1e3); + + MBinning binsthetatru("BinningThetatru"); + binsthetatru.SetEdges(7, -2.5, 32.5); + SetBinning((TH1*)aeff, &binsetru, &binsthetatru); + + const Int_t netru = aeff->GetNbinsX(); + const Int_t ntheta = aeff->GetNbinsY(); + + for (int j=1; j<=netru; j++) + { + for (int k=1; k<=ntheta; k++) + { + Double_t cont = 10000.0;; + ((TH1*)aeff)->SetBinContent(j,k,cont); + + Double_t dcont = 100.0; + ((TH1*)aeff)->SetBinError(j,k,dcont); + } + } + // *fLog << "Dummy aeff : netru =" << netru << ", ntheta = " << ntheta << endl; + //.................................... + + // number of Eunf and Var bins (histograms : fHUnfold, fHFlux) + const Int_t nEunf = fHFlux.GetNbinsX(); + const Int_t nVar = fHFlux.GetNbinsY(); + + // number of Etru and Theta bins (histogram *aeff of collection area) + + const Int_t nEtru = aeff->GetNbinsX(); + const Int_t nTheta = aeff->GetNbinsY(); + + //*fLog << "nEunf =" << nEunf << ", nVar = " << nVar << endl; + //*fLog << "nEtru =" << nEtru << ", nTheta = " << nTheta << endl; + + //................................................................... + // calculate effective collection area + // for the Eunf and Var bins of the histogram fHUnfold + // from the histogram *aeff, which has bins in Etru and Theta + // the result is the histogram fHAeff + // + + TH2D fHAeff; + fHAeff.Sumw2(); + SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold); + + Int_t errflag; + Double_t c0, c1, c2; + Double_t t1, t2, t3; + Double_t a1, a2, a3; + + Double_t *aeffbar; + aeffbar = new Double_t[nEtru]; + Double_t *daeffbar; + daeffbar = new Double_t[nEtru]; + + Double_t aeffEunfVar; + Double_t daeffEunfVar; + + + //------ start n loop ------ + for (int n=1; n<=nVar; n++) + { + Double_t Thetabar = thetabar->GetBinContent(n); + Double_t cosThetabar = cos(Thetabar); + + // determine Theta bins (k1, k2, k3) for interpolation in Theta + // k0 denotes the Theta bin from whicvh the error is copied + Int_t k0=0, k1=0, k2=0, k3=0; + for (int k=3; k<=nTheta; k++) + { + Double_t Thetalow = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(k); + if (Thetabar < Thetalow) + { + k1 = k-2; + k2 = k-1; + k3 = k; + k0 = k2; + break; + } + } + + if (k3 == 0) + { + k1 = nTheta-2; + k2 = nTheta-1; + k3 = nTheta; + k0 = k2; + } + + if (Thetabar < ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(2)) + k0 = 1; + else if (Thetabar > ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta)) + k0 = nTheta; + + Double_t Thetamin = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(1); + Double_t Thetamax = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta+1); + if (Thetabar < Thetamin || Thetabar > Thetamax) + { + *fLog << "MHFlux.cc : extrapolation in Theta; Thetabar = " << Thetabar + << ", Thetamin =" << Thetamin + << ", Thetamax =" << Thetamax << endl; + } + + //*fLog << "Var bin " << n << ":" << endl; + //*fLog << "Thetabar= " << Thetabar + // << ", k0= " << k0 + // << ", k1= " << k1 + // << ", k2= " << k2 + // << ", k3= " << k3 << endl; + + + // calculate effective collection area at Theta = Thetabar + // by quadratic interpolation in cos(Theta); + // do this for each bin of Etru + // + + for (int j=1; j<=nEtru; j++) + { + c0 = 0.0; + c1 = 0.0; + c2 = 0.0; + + t1 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k1) ); + t2 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k2) ); + t3 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k3) ); + + a1 = aeff->GetBinContent(j,k1); + a2 = aeff->GetBinContent(j,k2); + a3 = aeff->GetBinContent(j,k3); + + Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag); + aeffbar[j] = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar; + daeffbar[j] = aeff->GetBinError(j,k0); + + //*fLog << "Etru bin " << j << ": tbar= " << Thetabar + // << ", abar= " << aeffbar[j] + // << ", dabar= " << daeffbar[j] << endl; + } + + //--- start m loop --- + // calculate effective collection area at (E = Ebar, Theta = Thetabar) + // by quadratic interpolation in log10(Etru) + // do this for each bin of Eunf + // + for (int m=1; m<=nEunf; m++) + { + Double_t log10Ebar = 0.5 * + ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m) ) + +log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) ); + Double_t Ebar = pow(10.0, log10Ebar); + + // determine Etru bins (j1, j2, j3) for interpolation in E + // j0 denotes the Etru bin from which the error is copied + Int_t j0=0, j1=0, j2=0, j3=0; + + for (int j=3; j<=nEtru; j++) + { + Double_t Elow = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(j); + if (Ebar < Elow) + { + j1 = j-2; + j2 = j-1; + j3 = j; + j0 = j2; + break; + } + } + + if (j3 == 0) + { + j1 = nEtru-2; + j2 = nEtru-1; + j3 = nEtru; + j0 = j2; + } + + if (Ebar < ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(2)) + j0 = 1; + else if (Ebar > ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru)) + j0 = nEtru; + + Double_t Etrumin = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(1); + Double_t Etrumax = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru+1); + if (Ebar < Etrumin || Ebar > Etrumax) + { + *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar + << ", Etrumin =" << Etrumin + << ", Etrumax =" << Etrumax << endl; + } + + //*fLog << "Var bin " << n << ":" << endl; + //*fLog << "Ebar= " << Ebar + // << ", j1= " << j1 + // << ", j2= " << j2 + // << ", j3= " << j3 << endl; + + + c0=0.0; + c1=0.0; + c2=0.0; + + t1 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1) ) + +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1+1)) ); + + t2 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2) ) + +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2+1)) ); + + t3 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3) ) + +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3+1)) ); + + + a1 = aeffbar[j1]; + a2 = aeffbar[j2]; + a3 = aeffbar[j3]; + + Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag); + aeffEunfVar = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar); + daeffEunfVar = daeffbar[j0]; + + //*fLog << "Eunf bin " << m << ": Ebar= " << Ebar + // << ", aeffEunfVar = " << aeffEunfVar + // << ", daeffEunfVar = " << daeffEunfVar << endl; + + fHAeff.SetBinContent(m,n,aeffEunfVar); + fHAeff.SetBinError(m,n,daeffEunfVar); + } + //--- end m loop --- + } + //------ end n loop ------ + delete aeffbar; + + //................................................................... + // now calculate the absolute gamma flux + // + for (int m=1; m<=nEunf; m++) + { + Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m); + + for (int n=1; n<=nVar; n++) + { + Double_t Ngam = fHUnfold.GetBinContent(m,n); + Double_t dNgam = fHUnfold.GetBinError(m,n); + + Double_t Aeff = fHAeff.GetBinContent(m,n); + Double_t dAeff = fHAeff.GetBinError(m,n); + + Double_t Effon = teff->GetBinContent(n); + Double_t dEffon = teff->GetBinError(n); + + Double_t Cont, dCont; + if (Ngam > 0.0 && DeltaE > 0.0 && Effon > 0.0 && Aeff > 0.0) + { + Cont = Ngam / (DeltaE * Effon * Aeff); + dCont = Cont * sqrt( dNgam*dNgam / (Ngam*Ngam) + + dEffon*dEffon / (Effon*Effon) + + dAeff*dAeff / (Aeff*Aeff) ); + } + else + { + Cont = 1.e-20; + dCont = 1.e-20; + } + + fHFlux.SetBinContent(m,n,Cont); + fHFlux.SetBinError(m,n,dCont); + + //*fLog << "Eunf bin " << m << ", Var bin " << n + // << ": Ngam = " << Ngam << ", Flux = " + // << Cont << ", dFlux = " << dCont << endl; + //*fLog << ", DeltaE = " << DeltaE << ", Effon = " << Effon + // << ", Aeff = " << Aeff << endl; + } + } + + //.............................................. + // draw the differential photon flux vs. E-unf + // for the individual bins of the variable Var + + if (Draw == kTRUE) + { + for (int n=1; n<=nVar; n++) + { + TString strg0("Flux-"); + strg0 += fVarname; + TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E"); + + char txt[100]; + TString strg1("Photon flux vs. E-unfold for "); + TString strg2 = fVarname; + strg2 += "-bin %d"; + sprintf(txt, strg2, n); + TString strg3 = strg1 + txt; + + new TCanvas(txt, strg3); + // TCanvas &c = *MakeDefCanvas(txt, txt); + // gROOT->SetSelectedPad(NULL); + + gPad->SetLogx(); + + h.SetName(txt); + h.SetTitle(strg3); + h.SetXTitle("E-unfold [GeV] "); + h.SetYTitle("photons / (s m2 GeV)"); + h.DrawCopy(); + + // c.Modified(); + // c.Update(); + } + } + //........................ +} + +// ------------------------------------------------------------------------- +// +// Draw copies of the histograms +// +TObject *MHFlux::DrawClone(Option_t *opt) const +{ + TCanvas &c = *MakeDefCanvas("flux", "Orig - Unfold - Flux plots"); + c.Divide(2, 2); + + gROOT->SetSelectedPad(NULL); + + c.cd(1); + ((TH2*)&fHOrig)->DrawCopy(""); + gPad->SetLogx(); + + c.cd(2); + ((TH2*)&fHUnfold)->DrawCopy(""); + gPad->SetLogx(); + + c.cd(3); + ((TH2*)&fHFlux)->DrawCopy(""); + gPad->SetLogx(); + + c.Modified(); + c.Update(); + + return &c; +} + +// ------------------------------------------------------------------------- +// +// Draw the histograms +// +void MHFlux::Draw(Option_t *opt) +{ + if (!gPad) + MakeDefCanvas("flux", "orig-unfold-flux plots"); + + gPad->Divide(2,2); + + gPad->cd(1); + fHOrig.Draw(opt); + + gPad->cd(2); + fHUnfold.Draw(opt); + + gPad->cd(3); + fHFlux.Draw(opt); + + gPad->Modified(); + gPad->Update(); +} + +// ------------------------------------------------------------------------- +// +// Quadratic interpolation +// +// *** calculate the parameters of a parabula +// y = a + b*x + c*x**2 = F(x) +// such that yi = F(xi) for (i=1,3) +// +void MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3, + Double_t y1, Double_t y2, Double_t y3, + Double_t *a, Double_t *b, Double_t *c, Int_t *errflag) +{ + double ai11,ai12,ai13,ai21,ai22,ai23,ai31,ai32,ai33; + double det,det1; + //double yt1,yt2,yt3; + + det = x2*x3*x3 + x1*x2*x2 + x3*x1*x1 + - x2*x1*x1 - x3*x2*x2 - x1*x3*x3; + + if (det != 0.0) + { + *errflag = 0; + det1 = 1.0/det; + + ai11 = x2*x3*x3 - x3*x2*x2; + ai12 = x3*x1*x1 - x1*x3*x3; + ai13 = x1*x2*x2 - x2*x1*x1; + + ai21 = x2*x2 - x3*x3; + ai22 = x3*x3 - x1*x1; + ai23 = x1*x1 - x2*x2; + + ai31 = x3 - x2; + ai32 = x1 - x3; + ai33 = x2 - x1; + + *a = (ai11*y1 + ai12*y2 + ai13*y3) * det1; + *b = (ai21*y1 + ai22*y2 + ai23*y3) * det1; + *c = (ai31*y1 + ai32*y2 + ai33*y3) * det1; + + //yt1 = *a + *b * x1 + *c * x1*x1; + //yt2 = *a + *b * x2 + *c * x2*x2; + //yt3 = *a + *b * x3 + *c * x3*x3; + + //*fLog << "x1 = " << x1 << ", x2 = " << x2 << ", x3 = " << x3 << endl; + //*fLog << "y1 = " << y1 << ", y2 = " << y2 << ", y3 = " << y3 << endl; + //*fLog << "yt1 = " << yt1 << ", yt2 = " << yt2 + // << ", yt3 = " << yt3 << endl; + //*fLog << "*a = " << *a << ", *b = " << *b << ", *c= " << *c + // << ", *errflag = " << *errflag << endl; + + return; + } + + *errflag = 1; + *a = 0.0; + *b = 0.0; + *c = 0.0; + return; +} Index: trunk/MagicSoft/Mars/mhist/MHFlux.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFlux.h (revision 1415) +++ trunk/MagicSoft/Mars/mhist/MHFlux.h (revision 1415) @@ -0,0 +1,72 @@ +#ifndef MARS_MHFlux +#define MARS_MHFlux + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TH1 +#include "TH1.h" +#endif + +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +class TH2D; + +// base class MH is used because it defines "MakeDefCanvas" +// if base class MH is used one has to define the member function Fill +// because in MH Fill is set to zero +class MHFlux : public MH +{ +private: + const char *fVarname; + const char *fUnit; + + TH2D fHOrig; // original distribution of E-est + TH2D fHUnfold; // unfolded distribution of E-unfold + TH2D fHFlux; // absolute differential photon flux versus E-unfold + // all these plots for different bins of the variable (Theta or time) + +public: + MHFlux(const TH2D &h2d, const Bool_t Draw, + const char *varname=NULL, const char *unit=NULL); + + Bool_t Fill(const MParContainer *par); + + void Unfold(const Bool_t Draw); + void CalcFlux(const TH1D *teff, const TProfile *thetabar, + const TH2D *aeff, const Bool_t Draw); + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option="") const; + + const TH2D *GetHOrig() { return &fHOrig; } + const TH2D *GetHUnfold() { return &fHUnfold; } + const TH2D *GetHFlux() { return &fHFlux; } + + void Parab(double x1, double x2, double x3, + double y1, double y2, double y3, + double *a, double *b, double *c, int *errflag); + + ClassDef(MHFlux, 1) //2D-plots (original, unfolded, flux) +}; + +#endif + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhist/Makefile (revision 1414) +++ trunk/MagicSoft/Mars/mhist/Makefile (revision 1415) @@ -44,6 +44,5 @@ MHAlphaEnergyTime.cc \ MHAlphaEnergyTheta.cc \ - MHEffOnTimeTime.cc \ - MHEffOnTimeTheta.cc \ + MHEffOnTime.cc \ MHTimeDiffTime.cc \ MHTimeDiffTheta.cc \ @@ -55,10 +54,11 @@ MHMcEfficiencyEnergy.cc \ MHMcEnergyImpact.cc \ + MHMcRate.cc \ + MHThetabarTime.cc \ + MHThetabarTheta.cc \ MHMcEnergyMigration.cc \ - MHThetabarTime.cc \ - MHMcRate.cc \ - MHMcIntRate.cc \ - MHMcDifRate.cc - + MHGamma.cc \ + MHFlux.cc + SRCS = $(SRCFILES) @@ -77,2 +77,12 @@ # @endcode + + + + + + + + + +