/* ======================================================================== *\ ! ! * ! * 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"); }