Index: branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C =================================================================== --- branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C (revision 18094) +++ branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C (revision 18100) @@ -23,5 +23,6 @@ write->AddContainer("MTime", "Events", kFALSE); write->AddContainer("MMcEvt", "Events", kFALSE); - write->AddContainer("MEnergyEst", "Events", kFALSE); + //write->AddContainer("MEnergyEst", "Events", kFALSE); + //write->AddContainer("MEnergyTable", "Events", kFALSE); write->AddContainer("DataType", "Events"); write->AddContainer("Weight", "Events"); @@ -170,4 +171,8 @@ energyest.SetRule("MHillas.fSize^0.815*pow(10, 0.0199*(MHillasSrc.fDist*1.17e-2 - (-6.95))^2*(MHillasSrc.fDist*1.17e-2>(-6.95)) + 0.5997 *(MHillasSrc.fDist*1.17e-2 - 0.7568)^2 * (MHillasSrc.fDist*1.17e-2>0.7568) )"); // energyest.SetRule("1737 * (MHillas.fSize / 300.)^0.810"); + + //MEnergyTable energytable; + // Tables available from GH + //energytable.SetTableFile("Table_standard_test.root"); // Cuts before writing @@ -283,4 +288,5 @@ tlist2.AddToList(&cuts); tlist2.AddToList(&energyest); + //tlist2.AddToList(&energytable); tlist2.AddToList(write0); tlist2.AddToList(&cont0); // Post/Pre shown cuts Index: branches/Mars_McMismatchStudy/manalysis/AnalysisLinkDef.h =================================================================== --- branches/Mars_McMismatchStudy/manalysis/AnalysisLinkDef.h (revision 18094) +++ branches/Mars_McMismatchStudy/manalysis/AnalysisLinkDef.h (revision 18100) @@ -25,3 +25,5 @@ #pragma link C++ class MMcCalibrationUpdate+; +#pragma link C++ class MEnergyTable+; + #endif Index: branches/Mars_McMismatchStudy/manalysis/MEnergyTable.cc =================================================================== --- branches/Mars_McMismatchStudy/manalysis/MEnergyTable.cc (revision 18100) +++ branches/Mars_McMismatchStudy/manalysis/MEnergyTable.cc (revision 18100) @@ -0,0 +1,207 @@ +/* ======================================================================== *\ +! +! * +! * 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 9/2002 +! Wolfgang Wittek 7/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MEnergyTable // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MEnergyTable.h" + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MMcEvt.hxx" +#include "MPointingPos.h" +#include "MHMatrix.h" +#include "MHillas.h" +#include "MHillasSrc.h" +#include "MNewImagePar.h" +#include "MEnergyEstimate.h" +#include "MParameters.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "TFile.h" + +ClassImp(MEnergyTable); + +// -------------------------------------------------------------------------- +// +// Default constructor. Give the name of the parameter container (MHillas) +// storing width, length and size (Default="MHillas"). +// For the Zenith Angle MMcEvt.fTheta is used. +// +MEnergyTable::MEnergyTable(const char *name, const char *title) +{ + + fName = name ? name : "MEnergyTable"; + fTitle = title ? title : "Task to estimate the energy using MC Tables"; + + fNameHillasSrc = "MHillasSrc"; + fNameHillas = "MHillas"; + fNameMC = "MMcEvt"; + fNamePoint = "MPointingPos"; + //fNameParameter = "MParameterD"; + fNameParameter = "MEnergyTable"; + + SetNameParameter("MEnergyTable"); +} + + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MEnergyTable::~MEnergyTable() +{ + + delete fHillasSrc; + delete fHillas; + delete fParameter; + delete fMC; + delete fPointing; + + fTableFile->Close(); + +} + +Int_t MEnergyTable::PreProcess(MParList *pList) +{ + + fHillasSrc = (MHillasSrc*)pList->FindCreateObj("MHillasSrc", AddSerialNumber(fNameHillasSrc)); + if (!fHillasSrc) + return kFALSE; + + fHillas = (MHillas*)pList->FindCreateObj("MHillas", AddSerialNumber(fNameHillas)); + if (!fHillas) + return kFALSE; + + fPointing = (MPointingPos*)pList->FindCreateObj("MPointingPos", AddSerialNumber(fNamePoint)); + if (!fPointing) + return kFALSE; + + fMC = (MMcEvt*)pList->FindCreateObj("MMcEvt", AddSerialNumber(fNameMC)); + if (!fMC) + return kFALSE; + + fParameter = (MParameterD*)pList->FindCreateObj("MParameterD", fNameParameter); + if (!fParameter) + return kFALSE; + + + return kTRUE; + +} + +Int_t MEnergyTable::Process() +{ + + Double_t fDist = fHillasSrc->GetDist(); + Double_t fSize = fHillas->GetSize(); + fZenith = fPointing->GetZd(); + fThreshold = 305.; //TODO + + iZenith = 0; + iThreshold = 0; + + for( int i = 0; i < 2; i++ ) + if( fZenith >= iZen[i] && fZenith < iZen[i+1] ) iZenith = i; + + for( int i = 0; i < 2; i++ ) + if( fThreshold >= iThr[i] && fThreshold < iThr[i+1] ) iThreshold = i; + + if( iZenith != -999 && iThreshold != -999 ) + { + +// Cannot get Threshold Yet so doing 1D interpolation on Zenith angle +// Between i and i+1 bin + + iX = hEnergy[iZenith][0]->GetXaxis()->FindBin(log10(fSize)); + iY = hEnergy[iZenith][0]->GetYaxis()->FindBin(fDist); + fEnergy1 = hEnergy[iZenith][0]->GetBinContent(iX,iY); + + iX = hEnergy[iZenith+1][0]->GetXaxis()->FindBin(log10(fSize)); + iY = hEnergy[iZenith+1][0]->GetYaxis()->FindBin(fDist); + fEnergy2 = hEnergy[iZenith+1][0]->GetBinContent(iX,iY); + + fEnergyRecon = fEnergy2-fEnergy1; + fEnergyRecon /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]); + fEnergyRecon *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]); + fEnergyRecon += fEnergy1; + + } + + fParameter->SetVal(fEnergyRecon); + fParameter->SetReadyToSave(); + return kTRUE; + +} + +Int_t MEnergyTable::SetTableFile( const TString sTableFile ) +{ + + + //TFile *fTableFile = new TFile(sTableFile.c_str(),"OPEN"); + fTableFile = new TFile(sTableFile,"OPEN"); + if( fTableFile->IsZombie() ) + { + cout << "Cannot open Table File: " << sTableFile << endl; + return kFALSE; + } + + iZen[0] = 5; + iZen[1] = 30; + iZen[2] = 45; + iThr[0] = 300; + iThr[1] = 405; + iThr[2] = 500; + + for( int i = 0; i < 3; i++ ) + { + for( int j = 0; j < 3; j++ ) + { + sprintf(sHistDir,"z%.2d/t%.2d",iZen[i],iThr[j]); + sprintf(sHistName,"z%.2d/t%.2d/hEnergyAverage_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]); + if( fTableFile->GetDirectory(sHistDir) ) + { + fTableFile->cd(sHistDir); + cout << "Adding: " << i << " " << j << " " << sHistDir << endl; + htemp.push_back( (TH2F*)fTableFile->Get(sHistName) ); + } else + { + cout << "Error Table File not complete!" << endl; + return kFALSE; + } + } + hEnergy.push_back(htemp); + htemp.clear(); + } + + return kTRUE; + +} Index: branches/Mars_McMismatchStudy/manalysis/MEnergyTable.h =================================================================== --- branches/Mars_McMismatchStudy/manalysis/MEnergyTable.h (revision 18100) +++ branches/Mars_McMismatchStudy/manalysis/MEnergyTable.h (revision 18100) @@ -0,0 +1,95 @@ +#ifndef MARS_MEnergyTable +#define MARS_MEnergyTable + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +#include +#include "TH2F.h" + +#define D2R 3.1415/180. + +#ifndef MARS_MParameterCalc +#include "MParameterCalc.h" +#endif + +using namespace std; + +class MHillasSrc; +class MEnergyEst; +class MHillas; +class MNewImagePar; +class MMcEvt; +class MHMatrix; +class MPointingPos; +class MParameterD; + +class MEnergyTable : public MParameterCalc +// : public MTask +{ +private: + + MParameterD *fParameter; + MHillasSrc *fHillasSrc; + MHillas *fHillas; + MPointingPos *fPointing; + MMcEvt *fMC; + + TString fNameHillasSrc; + TString fNameHillas; + TString fNamePoint; + TString fNameMC; + TString fNameParameter; + + Double_t fMm2Deg; + + Int_t PreProcess(MParList *plist); + Int_t Process(); + + TFile *fTableFile; + + char sHistDir[500]; + char sHistName[500]; + int iZen[3]; + int iThr[3]; + + vector htemp; + vector< vector > hEnergy; + + Float_t fZenith; + float fThreshold; + Int_t iZenith; + Int_t iThreshold; + Int_t iX; + Int_t iY; + float fEnergy1; + float fEnergy2; + float fEnergyRecon; + +public: + + // TODO constructer will need some options? + //MEnergyTable(const char *hil="MHillas", const char *name=NULL, const char *title=NULL); + MEnergyTable(const char *name=NULL, const char *title=NULL); + ~MEnergyTable(); + + Int_t SetTableFile( const TString sTableFile ); + + ClassDef(MEnergyTable, 0) // Task to estimate the energy + +}; + +#endif + + + + + + + +