Index: trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyThetaONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyThetaONOFF.cc (revision 6967) +++ trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyThetaONOFF.cc (revision 6967) @@ -0,0 +1,420 @@ +/* ======================================================================== *\ +! +! * +! * 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 12/2004 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHExcessEnergyThetaONOFF // +// // +// 3D-histogram in alpha vs. E-est and Theta // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHExcessEnergyThetaONOFF.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" + +#include + +#include "MHMyFindSignificanceONOFF.h" + + +ClassImp(MHExcessEnergyThetaONOFF); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHExcessEnergyThetaONOFF::MHExcessEnergyThetaONOFF(const char *name, const char *title) + : fHist("","",10,0,1, 10,0,1) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHExcessEnergyThetaONOFF"; + 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 MHExcessEnergyThetaONOFF::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 MHExcessEnergyThetaONOFF::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 MHExcessEnergyThetaONOFF::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 MHExcessEnergyThetaONOFF::Calc(MHAlphaEnergyTheta* hAlphaON, MHAlphaEnergyTheta* hAlphaOFF) +{ + *fLog << endl; + fLog->Separator("Excess Events Calculation"); + + + TH3D* aetON = (TH3D*)hAlphaON->GetHist(); + TH3D* aetOFF = (TH3D*)hAlphaOFF->GetHist(); + + const TAxis* axisEnergy = aetON->GetXaxis(); + const TAxis* axisTheta = aetON->GetYaxis(); + const Int_t energyBins = aetON->GetXaxis()->GetNbins(); + const Int_t thetaBins = aetON->GetYaxis()->GetNbins();; + + MH::SetBinning(&fHist, axisEnergy, axisTheta); + + + + // + // Loop over Theta bins + // + for(Int_t iy=1; iy<=thetaBins; iy++) + { + TCanvas* c= new TCanvas(Form("zbin %d",iy),Form("zbin %d",iy)); + c->Divide(4,3); + + // + // Loop over Energy bins + // + for(Int_t ix=1; ix<=energyBins; ix++) + { + const Double_t e = aetON->GetXaxis()->GetBinCenter(ix); + const Double_t e1 = aetON->GetXaxis()->GetBinLowEdge(ix); + const Double_t e2 = aetON->GetXaxis()->GetBinLowEdge(ix+1); + + // + // Get Alpha plot for that Theta & Energy bin + // + TH1* halphaON = aetON->ProjectionZ("AlphaTempON",ix,ix,iy,iy); + halphaON->SetTitle(Form("Energy Bin = (%.0f, %.0f) GeV",e1,e2)); + halphaON->SetDirectory(NULL); + + TH1* halphaOFF = aetOFF->ProjectionZ("AlphaTempOFF",ix,ix,iy,iy); + halphaOFF->SetTitle(Form("Energy Bin = (%.0f, %.0f) GeV",e1,e2)); + halphaOFF->SetDirectory(NULL); + + TCanvas *c = new TCanvas; + halphaOFF->Draw(); + + + // + // Fit Alpha plot to get Nex + // + MHMyFindSignificanceONOFF findsig; + findsig.SetRebin(kFALSE); + //findsig.SetRebin(kTRUE); + double alphasig=25; + + // 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 = 10; +// break; +// case 7: +// alphasig = 10; +// break; +// case 8: +// alphasig = 10; +// break; +// case 9: +// alphasig = 10; +// break; +// case 10: + +// alphasig = 10; +// break; +// } + + const Double_t AlphaMinON = 30; + const Double_t AlphaMaxON = 90; + const Double_t AlphaMinOFF = 0; + const Double_t AlphaMaxOFF = AlphaMaxON; + const Double_t Degree = 4; + + fLog->SetNullOutput(kTRUE); + + + Bool_t rc = findsig.FindSigmaONOFF(halphaON, halphaOFF, alphasig, AlphaMinON,AlphaMaxON,Degree); + + // Double_t SigmaGauss = findsig.GetSigmaGauss(); + + // alphasig = SigmaGauss*3; +// rc = findsig.FindSigma(halpha, 25,90, degree, alphasig, 0,1,0); + + fLog->SetNullOutput(kFALSE); + + + + if(ix<=20) + { + c->cd(ix); + gPad->SetGridx(); + gPad->SetGridy(); + + findsig.Draw("all"); + } + + 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"); + + TCanvas *c2 = new TCanvas; + c2->SetLogx(); + + TH1* h= fHist.ProjectionX(); + h->SetStats(0); + h->SetTitle("No. Excess events Vs. Energy"); + h->SetXTitle("E [GeV]"); + h->SetYTitle("# Excess events"); + + h->Draw("e1"); + + + ///////////// + TCanvas *c3 = new TCanvas; + + + Double_t res[10] ={0.224281,0.30083, 0.397355, 0.455888, 0.482985, 0.479416, 0.441429, 0.424654, 0.366658, 0.287384}; + + TGraphErrors* graphIntegral = new TGraphErrors(h); + + for(int i=0; iGetN(); i++) + { + Double_t errorX = graphIntegral->GetErrorX(i); + Double_t errorY = graphIntegral->GetErrorY(i); + errorX *= res[i]; + + graphIntegral->SetPointError(i,errorX,errorY); + } + + + c3->SetLogx(); + c3->SetLogy(); + c3->SetGridx(); + c3->SetGridy(); + + graphIntegral->SetMarkerStyle(8); + graphIntegral->GetHistogram() ->SetXTitle("E [GeV]"); + graphIntegral->GetHistogram() ->SetYTitle("# Excess events"); + graphIntegral->Draw("Ap"); + +} Index: trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyThetaONOFF.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyThetaONOFF.h (revision 6967) +++ trunk/MagicSoft/Mars/mtemp/mucm/classes/MHExcessEnergyThetaONOFF.h (revision 6967) @@ -0,0 +1,44 @@ +#ifndef MARS_MHExcessEnergyThetaONOFF +#define MARS_MHExcessEnergyThetaONOFF + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH3 +#include +#endif + +#ifndef MARS_MHExcessEnergyTheta +#include "MHExcessEnergyTheta.h" +#endif + +class MParList; +class MHAlphaEnergyTheta; + +class MHExcessEnergyThetaONOFF : public MHExcessEnergyTheta +{ + private: + + TH2D fHist; + + Bool_t SetupFill(const MParList *pList); + Bool_t ReInit(MParList *pList); + + + public: + + MHExcessEnergyThetaONOFF(const char *name=NULL, const char *title=NULL); + + void Calc(MHAlphaEnergyTheta* hAlphaON, MHAlphaEnergyTheta* hAlphaOFF); + + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + + TH2D *GetHist() { return &fHist; } + + void Draw(Option_t *option=""); + + ClassDef(MHExcessEnergyThetaONOFF, 1) //2D-histogram in number of excess events vs. Energy and theta +}; + +#endif