Index: /trunk/MagicSoft/Mars/mtemp/Changelog =================================================================== --- /trunk/MagicSoft/Mars/mtemp/Changelog (revision 5449) +++ /trunk/MagicSoft/Mars/mtemp/Changelog (revision 5450) @@ -18,4 +18,10 @@ -*-*- END OF LINE -*-*- + +2004/11/22: Marcos Lopez + * mtemp/mucm/classes + - add new classes MHFlux.{h,cc} and MHExcessEnergyTheta.{h,cc} for the + flux calculation. + 2004/10/05: Javier Rico Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyTheta.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyTheta.cc (revision 5450) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyTheta.cc (revision 5450) @@ -0,0 +1,326 @@ +/* ======================================================================== *\ +! +! * +! * 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): Marcos Lopez 11/2004 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHExcessEnergyTheta // +// // +// 3D-histogram in alpha vs. E-est and Theta // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHExcessEnergyTheta.h" + +#include +#include +#include +#include + +#include "MMcEvt.hxx" +#include "MHillasSrc.h" +#include "MEnergyEst.h" +#include "MPointingPos.h" +#include "MRawRunHeader.h" +#include "MHAlphaEnergyTheta.h" +#include "MHFindSignificance.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHExcessEnergyTheta); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHExcessEnergyTheta::MHExcessEnergyTheta(const char *name, const char *title) + : fHist("","",10,0,1, 10,0,1) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHExcessEnergyTheta"; + fTitle = title ? title : "# Excess events vs. E and Theta"; + + + fHist.SetDirectory(NULL); + + fHist.SetTitle("# Excess events vs. E and Theta"); + fHist.SetXTitle("E [GeV]"); + fHist.SetYTitle("\\Theta [\\circ]"); + fHist.SetZTitle("# Excess events"); +} + +// -------------------------------------------------------------------------- +// +// Set binnings and prepare filling of the histogram +// +Bool_t MHExcessEnergyTheta::SetupFill(const MParList *plist) +{ + + // fEnergy = (MEnergyEst*)plist->FindObject("MEnergyEst"); +// if (!fEnergy) +// { +// *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl; +// return kFALSE; +// } + + +// // +// // Binning +// // +// MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); +// MBinning* binsalphaflux = (MBinning*)plist->FindObject("BinningAlphaFlux"); +// MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); +// if (!binsenergy || !binsalphaflux || !binstheta) +// { +// *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; +// return kFALSE; +// } + +// // SetBinning(&fHist, binsalphaflux, binsenergy, binstheta); +// SetBinning(&fHist, binsenergy, binstheta, binsalphaflux); + +// fHist.Sumw2(); + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Look in the parlist for MMcEvt or MPointingPos depending on the run type. +// +Bool_t MHExcessEnergyTheta::ReInit(MParList *pList) +{ + +// // This must be done in ReInit because in PreProcess the +// // runheaders are not available +// const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); +// if (!runheader) +// { +// *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; +// *fLog << warn << dbginf << "Warning - Assuming data type file...searching for MPointingPos..." << endl; +// } + + +// if (runheader && runheader->IsMonteCarloRun()) +// { +// fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); +// if (!fMcEvt) +// { +// *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; +// return kFALSE; +// } +// } +// else +// { +// fPointingPos = (MPointingPos*)pList->FindObject("MPointingPos"); +// if (!fPointingPos) +// { +// *fLog << err << dbginf << "PointingPos not found... aborting." << endl; +// return kFALSE; +// } +// } + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// Fill the histogram +// +Bool_t MHExcessEnergyTheta::Fill(const MParContainer *par, const Stat_t w) +{ + // MHillasSrc &hil = *(MHillasSrc*)par; + +// // fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(), +// // fMcEvt->GetTelescopeTheta()*kRad2Deg, w); + + +// const Double_t Zd = (fMcEvt) ? fMcEvt->GetTelescopeTheta()*kRad2Deg : fPointingPos->GetZd() ; + +// fHist.Fill(fEnergy->GetEnergy(), Zd, hil.GetAlpha(), w); + +// // *fLog<< Zd << " " << fEnergy->GetEnergy() << " " << hil.GetAlpha() << endl; + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// Get the Alpha distribution for each bin in energy and theta, +// fit it with MHFindSignificance to get the number of Excess events +// and it error. +// +void MHExcessEnergyTheta::Calc(MHAlphaEnergyTheta* hAlpha) +{ + TH3D* aet = (TH3D*)hAlpha->GetHist(); + + const TAxis* axisEnergy = aet->GetXaxis(); + const TAxis* axisTheta = aet->GetYaxis(); + const Int_t energyBins = aet->GetXaxis()->GetNbins(); + const Int_t thetaBins = aet->GetYaxis()->GetNbins();; + + MH::SetBinning(&fHist, axisEnergy, axisTheta); + + for(Int_t iy=1; iy<=thetaBins; iy++) + { + TCanvas* c= new TCanvas(Form("zbin %d",iy),Form("zbin %d",iy)); + c->Divide(4,3); + + for(Int_t ix=1; ix<=energyBins; ix++) + { const Double_t e = aet->GetXaxis()->GetBinCenter(ix); + const Double_t e1 = aet->GetXaxis()->GetBinLowEdge(ix); + const Double_t e2 = aet->GetXaxis()->GetBinLowEdge(ix+1); + + //TH1* halpha = aet->ProjectionZ("AlphaTemp",ix,ix,iy,iy); + TH1* halpha = aet->ProjectionZ("AlphaTemp",ix,ix); + halpha->SetTitle(Form("Energy Bin = (%.0f, %.0f) GeV",e1,e2)); + halpha->SetDirectory(NULL); + + MHFindSignificance findsig; + //findsig.SetRebin(kFALSE); + double alphasig=20; + + if(ix<5) + alphasig = 20; + else + alphasig = 20; + + switch (ix) + { + case 1: + alphasig = 20; + break; + case 2: + alphasig = 20; + break; + case 3: + alphasig = 20; + break; + case 4: + alphasig = 20; + break; + case 5: + alphasig = 20; + break; + case 6: + alphasig = 20; + break; + case 7: + alphasig = 15; + break; + case 8: + alphasig = 10; + break; + case 9: + alphasig = 5; + break; + case 10: + alphasig = 5; + break; + } + + + fLog->SetNullOutput(kTRUE); + + Bool_t rc = findsig.FindSigma(halpha, 30,90, 2, alphasig, 0,1,0); + Double_t SigmaGauss = findsig.GetSigmaGauss(); + + // alphasig = SigmaGauss*3; + +// rc = findsig.FindSigma(halpha, 30,90, 2, alphasig, 0,1,0); + + fLog->SetNullOutput(kFALSE); + + if(ix<=20) + { + c->cd(ix); + findsig.DrawFit("brief"); + } + + Double_t Nex = 0; + Double_t dNex = 0; + Double_t Sig = 0; + + if(!rc) + { + cout << "ERROR: Fit not possible. " << endl; + } + // else + { + SigmaGauss = findsig.GetSigmaGauss(); + Nex = findsig.GetNex(); + dNex = findsig.GetdNex(); + Sig = findsig.GetSigLiMa(); + } + // Avoid to have a negative number of excess events, which + // will give a negative flux which doesn't make sense + if(Nex<0) + { + Nex = 0; + dNex = 0; + Sig = 0; + } + + fHist.SetBinContent(ix,iy, Nex); + fHist.SetBinError(ix,iy, dNex); + + + // Print + *fLog<< endl + << " Theta Bin = " << iy << endl + << " Energy bin width E = " << ix << endl + << " Alphasig = " <SetLogx(); + c1->SetLogz(); + + fHist.SetStats(0); + fHist.Draw("lego"); +} Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyTheta.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyTheta.h (revision 5450) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyTheta.h (revision 5450) @@ -0,0 +1,40 @@ +#ifndef MARS_MHExcessEnergyTheta +#define MARS_MHExcessEnergyTheta + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH3 +#include +#endif + + +class MParList; +class MHAlphaEnergyTheta; + +class MHExcessEnergyTheta : public MH +{ +private: + + TH2D fHist; + + Bool_t SetupFill(const MParList *pList); + Bool_t ReInit(MParList *pList); + + +public: + + MHExcessEnergyTheta(const char *name=NULL, const char *title=NULL); + + void Calc(MHAlphaEnergyTheta* hAlpha); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + + TH2D *GetHist() { return &fHist; } + + void Draw(Option_t *option=""); + + ClassDef(MHExcessEnergyTheta, 1) //2D-histogram in number of excess events vs. Energy and theta +}; + +#endif Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHFlux.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHFlux.cc (revision 5450) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHFlux.cc (revision 5450) @@ -0,0 +1,549 @@ +/* ======================================================================== *\ +! +! * +! * 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): Marcos Lopex 11/2004 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHFlux // +// // +// 3D-histogram in alpha vs. E-est and Theta // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHFlux.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "MHillasSrc.h" +#include "MEnergyEst.h" +#include "MPointingPos.h" +#include "MRawRunHeader.h" + +#include "MHExcessEnergyTheta.h" +#include "MHMcCollectionArea.h" +#include "MHEffectiveOnTime.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHFlux); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHFlux::MHFlux(const char *name, const char *title) + : fHist("","",10,0,1, 10,0,1), fAverageFlux("","",1,0,1) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHFlux"; + fTitle = title ? title : "Flux vs. E and Theta"; + + + fHist.SetDirectory(NULL); + fHist.SetName("Flux vs. E and Theta"); + fHist.SetTitle("Flux vs. E and Theta"); + fHist.SetXTitle("E [GeV]"); + fHist.SetYTitle("\\Theta [\\circ]"); + fHist.SetZTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]"); + + fAverageFlux.SetDirectory(NULL); + fAverageFlux.SetName("Average Flux"); + fAverageFlux.SetTitle("Average Flux"); + fAverageFlux.SetXTitle("E [GeV]"); + fAverageFlux.SetYTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]"); +} + +// -------------------------------------------------------------------------- +// +// Set binnings and prepare filling of the histogram +// +Bool_t MHFlux::SetupFill(const MParList *plist) +{ + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Look in the parlist for MMcEvt or MPointingPos depending on the run type. +// +Bool_t MHFlux::ReInit(MParList *pList) +{ + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// Fill the histogram +// +Bool_t MHFlux::Fill(const MParContainer *par, const Stat_t w) +{ + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Calc +// +void MHFlux::Calc(MHExcessEnergyTheta* hExcess, MHMcCollectionArea* hColArea, MHEffectiveOnTime* hEffTime) +{ + const TH2D* hex = hExcess->GetHist(); + const TH1D* hca = hColArea->GetHist(); + const TH1D* heot = &hEffTime->GetHEffOnTheta(); + + TAxis* axisEnergy = hex->GetXaxis(); + const TAxis* axisTheta = hex->GetYaxis(); + const Int_t energyBins = hex->GetXaxis()->GetNbins(); + const Int_t thetaBins = hex->GetYaxis()->GetNbins();; + + MH::SetBinning(&fHist,axisEnergy, axisTheta); + + + // + // Calculate flux for each energy and theta + // + for (Int_t iy=1; iy<=thetaBins; iy++) // loop on theta + { + // Get Effective Time [sec] and its error + Double_t t = heot->GetBinContent(iy); + const Double_t dt = heot->GetBinError(iy); + + if (t < 1e-3) + t = 0.0; + + for (Int_t ix=1; ix<=energyBins; ix++) + { + const Double_t n = hex->GetBinContent(ix,iy); + const Double_t dn = hex->GetBinError(ix,iy); + + // Get AreaEff and its error + Double_t energy = axisEnergy->GetBinCenter(ix); + Int_t bin = hca->GetXaxis()->FindBin(energy); + + // Get NumberExcessEvents and its error + const Double_t a = hca->GetBinContent(bin)*1e4; //cm^2 + const Double_t da = hca->GetBinError(bin) *1e4; //cm^2 + + // energy bin width in TeV + const Double_t en = axisEnergy->GetBinWidth(ix)*1e-3; //TeV + + // + // Check that we can calculate the flux for the current bin + // + if (t==0) + cout << "No_Ton "; + if (a==0) + cout << "No_Aeff "; + if (n==0) + cout << "No_Events "; + if ((t == 0) || (a == 0) || (n == 0)) { + cout << endl; + continue; + } + + // + // Flux calculation and its error + // + const Double_t flux = n/(en*t*a); + + // error propagation formula + const Double_t errN = dn/(en*a*t); + const Double_t errA = da * n/(en*t*a*a); + const Double_t errT = dt * n/(en*a*t*t); + const Double_t error = sqrt(errN*errN + errA*errA + errT*errT); + + cout << dn << " " << en << " " << a << " " << t << endl; + + fHist.SetBinContent(ix,iy,flux); + fHist.SetBinError(ix,iy,error); + + } //energy + } //theta + fHist.Print("all"); +} + + + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHFlux::Draw(Option_t *opt) +{ + + // -------------------------------------------------------------------- + // + // Draw lego plot of Flux vs. Energy and Theta + // + TCanvas *c1 = new TCanvas("Flux vs. E and Theta","Flux vs. E and Theta"); + c1->SetLogx(); + c1->SetLogz(); + fHist.SetStats(0); + fHist.Draw("lego"); + + + // -------------------------------------------------------------------- + // + // Draw the Flux for each Theta bin + // + TCanvas *c2 = new TCanvas("Fluxes for each Theta bin","Fluxes for each Theta bin"); + c2->SetLogx(); + c2->SetLogy(); + c2->SetGridx(); + c2->SetGridy(); + + + THStack* hs = new THStack("Fluxes for each Theta bin","Fluxes for each Theta bin"); + + TLegend * leg = new TLegend(0.73,0.65,0.89,0.89); + + TAxis* yaxis = fHist.GetYaxis(); + const Int_t nbiny = fHist.GetYaxis()->GetNbins(); + + for(Int_t iy=1; iy<=nbiny; iy++) + { + + TH1D* h1 = fHist.ProjectionX(Form("%d",iy),iy,iy,"e"); //<----- Option e is very important, otherwise the errors are not copied + + if(h1->GetEntries()==0) + continue; + + h1->SetLineColor(iy); + hs->Add(h1,"e1"); + leg->AddEntry(h1,Form("\\theta = %.0f",yaxis->GetBinCenter(iy)),"l"); + +// TCanvas *c = new TCanvas(); +// c->SetLogx(); +// c->SetLogy(); +// h1->DrawCopy("e"); +// h1->Print("all"); + } + + + + // -------------------------------------------------------------------- + // + // Calculate and Draw the Flux average on Theta bins + // + fAverageFlux.SetStats(0); + + MH::SetBinning(&fAverageFlux,fHist.GetXaxis()); + + for(int ix=1; ix<=fHist.GetXaxis()->GetNbins(); ix++) // energy + { + Double_t sumw = 0; + Double_t sumcontents = 0; + Double_t sumerrors = 0; + + for(int iy=1; iy<=fHist.GetYaxis()->GetNbins(); iy++) // theta + { + Double_t weight = fHist.GetYaxis()->GetBinWidth(iy); + sumw += weight; + sumcontents += fHist.GetBinContent(ix,iy)*weight; + sumerrors += fHist.GetBinError(ix,iy)*weight; + } + + fAverageFlux.SetBinContent(ix,sumcontents/sumw); + fAverageFlux.SetBinError(ix,sumerrors/sumw); + } + + // for(int ix=1; ix<=fHist.GetXaxis()->GetNbins(); ix++) // energy +// { +// Double_t sumw = 0; +// Double_t sumcontents = 0; + +// cout << " energy bin "<GetNbins(); iy++) // theta +// { + +// Double_t bincontent = fHist.GetBinContent(ix,iy); +// Double_t binerror = fHist.GetBinError(ix,iy); + +// if( bincontent == 0 || binerror == 0 ) +// continue; +// cout << binerror << endl; + +// Double_t weight = 1/(binerror*binerror); +// sumw += weight; +// sumcontents += bincontent*weight; + +// cout << " theta bin " << fHist.GetBinContent(ix,iy)<< " " <Add(&fAverageFlux,"pe1"); + leg->AddEntry(&fAverageFlux,"Average on Theta","l"); + + + c2->cd(); + hs->Draw("nostack"); + leg->Draw(); + + + TCanvas *c3 = new TCanvas("Average Flux","Average Flux"); + c3->SetLogx(); + c3->SetLogy(); + c3->SetGridx(); + c3->SetGridy(); + fAverageFlux.Draw(); + fAverageFlux.Print("all"); + + // + // Fix the Average Flux to a power law + // + TF1* fluxfit = new TF1("f1","[0]*pow(x,-[1])",90,1500); + fluxfit->SetParNames("f0","a"); + fluxfit->SetParameter(0,5.10986e-05); + fluxfit->SetParameter(1,2.4); + fluxfit->SetTitle("Flux fit"); + fluxfit->SetLineColor(27); + fluxfit->SetLineWidth(3); + + fAverageFlux.Fit("f1","R"); + + + // + // Draw the Crab spectrum measured by HEGRA between 500 GeV and 80 TeV + // + TF1* CrabFlux = new TF1("CrabFlux","[0]*pow(x/1000.,-[1])",350,2000); + CrabFlux->SetParameter(0,2.83e-11); + CrabFlux->SetParameter(1,2.62); + CrabFlux->SetLineStyle(2); + CrabFlux->SetLineColor(4); + CrabFlux->Draw("same"); + + // + // Draw formula + // + TPaveText* func = new TPaveText(0.16, 0.22, 0.67, 0.28,"NDC"); + func->AddText(Form("#frac{dF}{dE} = %.2e * E^{-%.2f} [#frac{ph}{cm^{2} s TeV}]",fluxfit->GetParameter(0),fluxfit->GetParameter(1))); + func->SetFillStyle(0); + func->SetBorderSize(0); + func->Draw(); + + +// // +// // Draw "Preliminary" +// // +// TPaveText* lab = new TPaveText(0.33, 0.83, 0.68, 0.89,"NDC"); +// lab->AddText("preliminary"); +// lab->SetTextColor(2); +// lab->SetFillStyle(0); +// lab->SetBorderSize(0); +// lab->Draw(); + + + + // --------------------------------------------------------------------- + // + // Integral flux + // + TH1D *hIntegral = (TH1D*)fAverageFlux.Clone(); + hIntegral->GetListOfFunctions()->Clear(); + + Int_t nbinsx = fAverageFlux.GetNbinsX(); + + for(int i=1; i<=nbinsx; i++) + { + cout <<"Integral Flux: Binwidth:" << hIntegral->GetBinWidth(i) << endl; + + hIntegral->SetBinContent(i,hIntegral->GetBinContent(i)*hIntegral->GetBinWidth(i)*1e-3); + hIntegral->SetBinError(i,hIntegral->GetBinError(i)*hIntegral->GetBinWidth(i)*1e-3); + } + + + for(int i=nbinsx-1; i>=1; i--) + { + Double_t integralsofar = hIntegral->GetBinContent(i+1); + Double_t current = hIntegral->GetBinContent(i); + + Double_t currentE = hIntegral->GetBinError(i); + Double_t Esofar = hIntegral->GetBinError(i+1); + + hIntegral->SetBinContent(i,(current+integralsofar)); + hIntegral->SetBinError(i,TMath::Sqrt(currentE*currentE+Esofar*Esofar)); + } + + + hIntegral->SetTitle("Integral Flux"); + hIntegral->SetXTitle("E [GeV]"); + hIntegral->SetYTitle("Integral Flux [s^{-1} cm^{-2}]"); + + + TCanvas *c20 = new TCanvas(); + c20->SetLogx(); + c20->SetLogy(); + c20->SetGridx(); + c20->SetGridy(); + + hIntegral->Draw(); + + + + + // -------------------------------------------------------------------- + // + // E^2 * Flux + // + TH1D *hEscaledFlux = (TH1D*)fAverageFlux.Clone(); + hEscaledFlux->GetListOfFunctions()->Clear(); + + nbinsx = hEscaledFlux->GetNbinsX(); + + for(int i=1; i<=nbinsx; i++) + { + + Double_t energy = hEscaledFlux->GetBinLowEdge(i)*1e-3; // TeV + Double_t Flux = hEscaledFlux->GetBinContent(i); + Double_t dFlux = hEscaledFlux->GetBinError(i); + + hEscaledFlux->SetBinContent(i,energy*energy*Flux); + hEscaledFlux->SetBinError(i,energy*energy*dFlux); + } + + TCanvas *c40 = new TCanvas(); + c40->SetLogx(); + c40->SetLogy(); + c40->SetGridx(); + c40->SetGridy(); + + hEscaledFlux->SetTitle("Escaled Flux"); + hEscaledFlux ->SetXTitle("E [GeV]"); + hEscaledFlux ->SetYTitle("E^{2}*Flux [TeV s^{-1} cm^{-2}]"); + + hEscaledFlux->Draw(); + + +// // ----------------------------------------------------------------- +// // +// // Graph move a 30 % +// // +// TCanvas *c4 = new TCanvas(); +// c4->SetLogx(); +// c4->SetLogy(); +// c4->SetGridx(); +// c4->SetGridy(); + +// Int_t nbins = fAverageFlux.GetNbinsX(); +// TArrayD x(nbins),y(nbins),dx(nbins),dy(nbins); + +// for(int i=1; i<=nbins; i++) +// { +// x[i-1] = fAverageFlux.GetXaxis()->GetBinCenter(i)*.7; +// y[i-1] = fAverageFlux.GetBinContent(i); +// dx[i-1] = fAverageFlux.GetXaxis()->GetBinWidth(i)*0.62; +// dy[i-1] = fAverageFlux.GetBinError(i); +// } + +// TGraphErrors* gr = new TGraphErrors(fAverageFlux.GetNbinsX(), x.GetArray(), y.GetArray(), dx.GetArray(), dy.GetArray()); + +// gr->SetMarkerStyle(8); +// gr->Draw("Ap"); +// gr->Print("all"); + +// // +// TF1* fluxfit2 = new TF1("f2","[0]*pow(x,-[1])",70,2500); +// fluxfit2->SetParNames("f0","a"); +// fluxfit2->SetParameter(0,5.10986e-05); +// fluxfit2->SetParameter(1,2.4); +// fluxfit2->SetTitle("Flux fit"); +// fluxfit2->SetLineColor(27); +// fluxfit2->SetLineWidth(3); + + +// gr->Fit("f2","R"); + +// gr->SetTitle(""); +// gr->SetMaximum(1e-3); +// gr->SetMinimum(1e-12); + +// TLegend* leg2 = new TLegend(0.67,0.72,0.89,0.89); + +// leg2->AddEntry(gr,"MAGIC Sept. 2004","p"); + +// gr->SetMarkerStyle(8); +// // gr->SetLineColor(6); + +// // // +// // // Draw the Crab spectrum measured by HEGRA between 500 GeV and 80 TeV +// // // +// // TF1* CrabFlux = new TF1("CrabFlux","[0]*pow(x/1000.,-[1])",350,2000); +// // CrabFlux->SetParameter(0,2.83e-11); +// // CrabFlux->SetParameter(1,2.62); +// // CrabFlux->SetLineStyle(2); +// // CrabFlux->SetLineColor(4); +// // CrabFlux->Draw("same"); + +// leg2->AddEntry(CrabFlux,"HEGRA ApJ 614","l"); +// leg2->Draw(); +// lab->Draw(); + +// TPaveText* func2 = new TPaveText(0.16, 0.22, 0.67, 0.28,"NDC"); +// func2->AddText(Form("#frac{dF}{dE} = %.2e * E^{-%.2f} [#frac{ph}{cm^{2} s TeV}]",fluxfit2->GetParameter(0),fluxfit2->GetParameter(1))); +// func2->SetFillStyle(0); +// func2->SetBorderSize(0); +// func2->Draw(); + + +// gr->GetHistogram()->SetXTitle("E [GeV]"); +// gr->GetHistogram()->SetYTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]"); + +// TH1F* h = gr->GetHistogram(); +// TCanvas *c12 = new TCanvas(); +// h->Draw(); +// cout << "Integral flux = "<< h->Integral("width") << endl; +} Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHFlux.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHFlux.h (revision 5450) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHFlux.h (revision 5450) @@ -0,0 +1,51 @@ +#ifndef MARS_MHFlux +#define MARS_MHFlux + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH3 +#include +#endif + +#ifndef ROOT_TH1 +#include +#endif + + +class MParList; + +class MHExcessEnergyTheta; +class MHMcCollectionArea; +class MHEffectiveOnTime; + + +class MHFlux : public MH +{ +private: + + TH2D fHist; + TH1D fAverageFlux; + + Bool_t SetupFill(const MParList *pList); + Bool_t ReInit(MParList *pList); + + +public: + + MHFlux(const char *name=NULL, const char *title=NULL); + + void Calc(MHExcessEnergyTheta* hex, MHMcCollectionArea* area, MHEffectiveOnTime* hEffTime); + + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + + TH2D *GetHist() { return &fHist; } + TH1D *GetAverageFlux() { return &fAverageFlux; } + + void Draw(Option_t *option=""); + + ClassDef(MHFlux, 1) //2D-histogram of Flux vs. Energy and theta +}; + +#endif