/* ======================================================================== *\ ! ! * ! * 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 // // Energy Tables generated at discrete Zenith and threshold // // We then interpolate between them to estimate the energy // // // ///////////////////////////////////////////////////////////////////////////// #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"; fNameParameterL = "ScaledLength"; fNameParameterW = "ScaledWidth"; SetNameParameter("MEnergyTable"); } // -------------------------------------------------------------------------- // // Destructor. // MEnergyTable::~MEnergyTable() { delete fHillasSrc; delete fHillas; delete fParameter; delete fParameterL; delete fParameterW; 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; fParameter = (MParameterD*)pList->FindCreateObj("MParameterD", fNameParameter); if (!fParameter) return kFALSE; fParameterW = (MParameterD*)pList->FindCreateObj("MParameterD", fNameParameterW); if (!fParameterW) return kFALSE; fParameterL = (MParameterD*)pList->FindCreateObj("MParameterD", fNameParameterL); if (!fParameterL) 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); fLength1 = hLength[iZenith][0]->GetBinContent(iX,iY); fLengthSigma1 = hLength[iZenith][0]->GetBinContent(iX,iY); fWidth1 = hLength[iZenith][0]->GetBinContent(iX,iY); fWidthSigma1 = hLength[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); fLength2 = hLength[iZenith+1][0]->GetBinContent(iX,iY); fLengthSigma2 = hLength[iZenith+1][0]->GetBinContent(iX,iY); fWidth2 = hLength[iZenith+1][0]->GetBinContent(iX,iY); fWidthSigma2 = hLength[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; fLength = fLength2-fLength1; fLength /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]); fLength *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]); fLength += fLength1; fLengthSigma = fLengthSigma2-fLengthSigma1; fLengthSigma /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]); fLengthSigma *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]); fLengthSigma += fLengthSigma1; fWidth = fWidth2-fWidth1; fWidth /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]); fWidth *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]); fWidth += fWidth1; fWidthSigma = fWidthSigma2-fWidthSigma1; fWidthSigma /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]); fWidthSigma *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]); fWidthSigma += fWidthSigma1; RSW = (fHillas->GetWidth()-fWidth)/fWidthSigma; RSL = (fHillas->GetLength()-fLength)/fLengthSigma; } fParameter->SetVal(fEnergyRecon); fParameter->SetReadyToSave(); fParameterW->SetVal(RSW); fParameterW->SetReadyToSave(); fParameterL->SetVal(RSL); fParameterL->SetReadyToSave(); return kTRUE; } Int_t MEnergyTable::SetTableFile( const TString sTableFile ) { fTableFile = new TFile(sTableFile,"OPEN"); if( fTableFile->IsZombie() ) { cout << "Cannot open Table File: " << sTableFile << endl; return kFALSE; } // TODO 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]); sprintf(sHistNameL,"z%.2d/t%.2d/hLengthAverage_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]); sprintf(sHistNameLS,"z%.2d/t%.2d/hLengthSigma_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]); sprintf(sHistNameW,"z%.2d/t%.2d/hWidthAverage_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]); sprintf(sHistNameWS,"z%.2d/t%.2d/hWidthSigma_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]); if( fTableFile->GetDirectory(sHistDir) ) { fTableFile->cd(sHistDir); cout << "Adding: " << i << " " << j << " " << sHistDir << endl; htempE.push_back( (TH2F*)fTableFile->Get(sHistName) ); htempL.push_back( (TH2F*)fTableFile->Get(sHistNameL) ); htempLS.push_back( (TH2F*)fTableFile->Get(sHistNameLS) ); htempW.push_back( (TH2F*)fTableFile->Get(sHistNameW) ); htempWS.push_back( (TH2F*)fTableFile->Get(sHistNameWS) ); } else { cout << "Error Table File not complete!" << endl; return kFALSE; } } hEnergy.push_back(htempE); hLength.push_back(htempL); hLengthSigma.push_back(htempLS); hWidth.push_back(htempW); hWidthSigma.push_back(htempWS); htempE.clear(); htempL.clear(); htempLS.clear(); htempW.clear(); htempWS.clear(); } return kTRUE; }