Index: /trunk/MagicSoft/Mars/macros/CT1EgyEst.C =================================================================== --- /trunk/MagicSoft/Mars/macros/CT1EgyEst.C (revision 2271) +++ /trunk/MagicSoft/Mars/macros/CT1EgyEst.C (revision 2272) @@ -33,4 +33,5 @@ #include "MHMatrix.h" #include "MEnergyEstParam.h" +#include "MEnergyEstParamDanielMkn421.h" #include "MMatrixLoop.h" #include "MChisqEval.h" @@ -138,5 +139,5 @@ Optimize.SetOutFile(paramout); - /* + MFCT1SelFinal filterhadrons; filterhadrons.SetHadronnessName(hadronnessName); @@ -144,5 +145,5 @@ filterhadrons.SetInverted(); Optimize.SetEventFilter(&filterhadrons); - */ + Optimize.SetNevents(howMany); @@ -205,9 +206,13 @@ // - MEnergyEstParam eest2(hilName); + //MEnergyEstParam eest2(hilName); + //eest2.Add(hilSrcName); + + //eest2.SetCoeffA(parA); + //eest2.SetCoeffB(parB); + + // estimate energy using Daniel's parameters + MEnergyEstParamDanielMkn421 eest2(hilName); eest2.Add(hilSrcName); - - eest2.SetCoeffA(parA); - eest2.SetCoeffB(parB); @@ -238,6 +243,5 @@ // Select events // - - /* + cout << "Select events with hadronness < " << maxhadronness << " and |alpha| < " << maxalpha << endl; @@ -247,30 +251,13 @@ MContinue cont(&hcut2); - */ - - parlist2.AddToList(&tlist2); - - //******************************** - // Entries in MTaskList - - tlist2.AddToList(&read2); - // tlist2.AddToList(&cont); - tlist2.AddToList(&eest2); - - // + + //----------------------------------------------- // Create Object MHMcEnergyMigration containing useful histograms, // and task MHMcEnergyMigration to fill them: // - MHMcEnergyMigration mighist; - - parlist2.AddToList(&mighist); - - MFillH migfill(&mighist, "MMcEvt"); - - tlist2.AddToList(&migfill); - - parlist2.AddToList(binsE); - parlist2.AddToList(binsTheta); + //MHMcEnergyMigration mighist; + + MFillH migfill("MHMcEnergyMigration", "MMcEvt"); MBinning BinningDE("BinningDE"); @@ -279,7 +266,24 @@ BinningDE.SetEdges(60, -1.2, 1.2); BinningImpact.SetEdges(50, 0., 400.); + + //******************************** + // Entries in MParList + + parlist2.AddToList(&tlist2); + //parlist2.AddToList(&mighist); + parlist2.AddToList(binsE); + parlist2.AddToList(binsTheta); parlist2.AddToList(&BinningDE); parlist2.AddToList(&BinningImpact); + //******************************** + // Entries in MTaskList + + tlist2.AddToList(&read2); + tlist2.AddToList(&cont); + tlist2.AddToList(&eest2); + tlist2.AddToList(&migfill); + + //----------------------------------------------- cout << "Event loop was setup" << endl; MProgressBar bar; @@ -291,4 +295,15 @@ if (!evtloop2.Eventloop()) return; + + tlist2.PrintStatistics(0, kTRUE); + + MHMcEnergyMigration &mighist = + *(MHMcEnergyMigration*)parlist2->FindObject("MHMcEnergyMigration"); + if (!&mighist) + { + cout << "CT1EgyEst.C : object 'MHMcEnergyMigration' not found ... aborting" + << endl; + return; + } TString paramout(outPath); Index: /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2271) +++ /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2272) @@ -30,4 +30,5 @@ #pragma link C++ class MEnergyEstimate+; #pragma link C++ class MEnergyEstParam+; +#pragma link C++ class MEnergyEstParamDanielMkn421+; #pragma link C++ class MMatrixLoop+; Index: /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2271) +++ /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2272) @@ -144,4 +144,125 @@ fAlphaUp[7] = 0.; //--------------------------------- +} + +// -------------------------------------------------------------------------- +// +// Set the parameter values from vector 'par' +// +// Attention : it is assumed that there are (9*ncutpar) values +// +void MCT1SupercutsCalc::SetParams(Double_t *par) +{ + UInt_t ncutpar = fLengthUp.GetSize(); + UInt_t k0 = 0; + + TArrayD lup(ncutpar, par + k0); + SetLengthUp(lup); + k0 += ncutpar; + + TArrayD wup(ncutpar, par + k0); + SetWidthUp(wup); + k0 += ncutpar; + + TArrayD dup(ncutpar, par + k0); + SetDistUp(dup); + k0 += ncutpar; + + TArrayD llo(ncutpar, par + k0); + SetLengthLo(llo); + k0 += ncutpar; + + TArrayD wlo(ncutpar, par + k0); + SetWidthLo(wlo); + k0 += ncutpar; + + TArrayD dlo(ncutpar, par + k0); + SetDistLo(dlo); + k0 += ncutpar; + + TArrayD aup(ncutpar, par + k0); + SetAsymUp(aup); + k0 += ncutpar; + + TArrayD alo(ncutpar, par + k0); + SetAsymLo(alo); + k0 += ncutpar; + + TArrayD alphaup(ncutpar, par + k0); + SetAlphaUp(alphaup); +} + +// -------------------------------------------------------------------------- +// +// Get the parameter values +// +// Attention : it is assumed that there are (9*ncutpar) values +// +void MCT1SupercutsCalc::GetParams(Double_t *par) +{ + UInt_t ncutpar = fLengthUp.GetSize(); + UInt_t k0 = 0; + + for (UInt_t j=0; jAddColumn("MHillas.fMeanY"); fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist"); + fMap[7] = fMatrix->AddColumn("MHillasSrc.fAlpha"); } @@ -332,15 +454,3 @@ return kTRUE; } - - - - - - - - - - - - - +//========================================================================== Index: /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2271) +++ /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2272) @@ -34,5 +34,5 @@ Double_t fMm2Deg; - Int_t fMap[7]; + Int_t fMap[8]; MHMatrix *fMatrix; @@ -70,4 +70,7 @@ const char *name=NULL, const char *title=NULL); + void SetParams(Double_t *par); + void GetParams(Double_t *par); + void SetHadronnessName(const TString name) { fHadronnessName = name; } TString GetHadronnessName() const { return fHadronnessName; } @@ -100,2 +103,10 @@ #endif + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MEnergyEstParamDanielMkn421.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParamDanielMkn421.cc (revision 2272) +++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParamDanielMkn421.cc (revision 2272) @@ -0,0 +1,401 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MEnergyEstParamDanielMkn421 // +// // +// Task to estimate the energy using a parametrization. // +// Make sure, that all source dependant containers are based on the same // +// set of basic hillas parameters // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MEnergyEstParamDanielMkn421.h" + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MMcEvt.hxx" +#include "MHMatrix.h" +#include "MHillasSrc.h" +#include "MNewImagePar.h" +#include "MEnergyEst.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MEnergyEstParamDanielMkn421); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Initialize the coefficients with Daniel's values +// for Mkn421 +// +// +void MEnergyEstParamDanielMkn421::InitCoefficients() +{ + // parameters for the impact parameter + fA.Set(4); + fA[0] = 12.0322; // + fA[1] = 8.29911; // + fA[2] = 73.0699; // + fA[3] = 0.311375; + + + // parameters for the energy + fB.Set(6); + fB[0] = 1.23138; // + fB[1] = 0.00276497; // + fB[2] = -0.0110667; // + fB[3] = 7.47538; // + fB[4] = -1.25619; // + fB[5] = 0.627699; // +} + +// -------------------------------------------------------------------------- +// +// 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. +// +MEnergyEstParamDanielMkn421::MEnergyEstParamDanielMkn421(const char *hillas, const char *name, const char *title) + : fMatrix(NULL) +{ + fName = name ? name : "MEnergyEstParamDanielMkn421"; + fTitle = title ? title : "Task to estimate the energy (Daniel Mkn421)"; + + fHillasName = hillas; + + fPairs = new TList; + fEnergy = new TList; + fHillasSrc = new TList; + + fPairs->SetOwner(); + + InitCoefficients(); + + AddToBranchList("MMcEvt.fTelescopeTheta"); + AddToBranchList(fHillasName+".fSize"); + AddToBranchList(fHillasName+".fWidth"); + AddToBranchList(fHillasName+".fLength"); + AddToBranchList("MNewImagePar.fLeakage1"); +} + + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MEnergyEstParamDanielMkn421::~MEnergyEstParamDanielMkn421() +{ + delete fPairs; + delete fEnergy; + delete fHillasSrc; +} + +// -------------------------------------------------------------------------- +// +// Check for all necessary parameter containers. +// +Int_t MEnergyEstParamDanielMkn421::PreProcess(MParList *plist) +{ + const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam"); + if (!geom) + { + *fLog << warn << dbginf << "No Camera Geometry available" << endl; + return kFALSE; + } + else + { + fMm2Deg = geom->GetConvMm2Deg(); + } + + + if (!fMatrix) + { + fHillas = (MHillas*)plist->FindObject(fHillasName, "MHillas"); + if (!fHillas) + { + *fLog << err << dbginf << fHillasName << " [MHillas] not found... aborting." << endl; + return kFALSE; + } + + fMc = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMc) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fNewImagePar = (MNewImagePar*)plist->FindObject("MNewImagePar"); + if (!fNewImagePar) + { + *fLog << err << dbginf << " object 'MNewImagePar' not found... aborting." << endl; + return kFALSE; + } + + } + + + TObject *obj = NULL; + TIter Next(fPairs); + + while ((obj=Next())) + { + TObject *o = plist->FindCreateObj(obj->GetTitle(), "MEnergyEst"); + if (!o) + return kFALSE; + + if (!fEnergy->FindObject(obj->GetTitle())) + fEnergy->Add(o); + + if (fMatrix) + continue; + + o = plist->FindObject(obj->GetName(), "MHillasSrc"); + if (!o) + { + *fLog << err << dbginf << obj->GetName() << " not found... aborting." << endl; + return kFALSE; + } + fHillasSrc->Add(o); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set the five coefficients for the estimation of the impact parameter. +// Argument must ba a TArrayD of size 5. +// +void MEnergyEstParamDanielMkn421::SetCoeffA(const TArrayD &arr) +{ + if (arr.GetSize() != fA.GetSize()) + { + *fLog << err << dbginf << "Error - Wrong number of coefficients!" << endl; + return; + } + + fA = arr; +} + +// -------------------------------------------------------------------------- +// +// Set the seven coefficients for the estimation of the energy. +// Argument must ba a TArrayD of size 7. +// +void MEnergyEstParamDanielMkn421::SetCoeffB(const TArrayD &arr) +{ + if (arr.GetSize() != fB.GetSize()) + { + *fLog << err << dbginf << "Error - Wrong number of coefficients!" << endl; + return; + } + + fB = arr; +} + +// -------------------------------------------------------------------------- +// +// Set the twelve coefficients for the estimation of impact and energy. +// Argument must ba a TArrayD of size 12. +// +void MEnergyEstParamDanielMkn421::SetCoeff(const TArrayD &arr) +{ + if (arr.GetSize() != fA.GetSize()+fB.GetSize()) + { + *fLog << err << dbginf << "Error - Wrong number of coefficients!" << endl; + return; + } + + fA = TArrayD(fA.GetSize(), arr.GetArray()); + fB = TArrayD(fB.GetSize(), arr.GetArray() + fA.GetSize()); +} + +// -------------------------------------------------------------------------- +// +// Returns the mapped value from the Matrix +// +Double_t MEnergyEstParamDanielMkn421::GetVal(Int_t i) const +{ + return (*fMatrix)[fMap[i]]; +} + +// -------------------------------------------------------------------------- +// +// You can use this function if you want to use a MHMatrix instead of the +// given containers. This function adds all necessary columns to the +// given matrix. Afterward you should fill the matrix with the corresponding +// data (eg from a file by using MHMatrix::Fill). If you now loop +// through the matrix (eg using MMatrixLoop) MEnergyEstParamDanielMkn421::Process +// will take the values from the matrix instead of the containers. +// +void MEnergyEstParamDanielMkn421::InitMapping(MHMatrix *mat) +{ + if (fMatrix) + return; + + fMatrix = mat; + + fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta"); + fMap[1] = fMatrix->AddColumn(fHillasName+".fWidth"); + fMap[2] = fMatrix->AddColumn(fHillasName+".fLength"); + fMap[3] = fMatrix->AddColumn(fHillasName+".fSize"); + fMap[4] = fMatrix->AddColumn("MNewImagePar.fLeakage1"); + + Int_t col = 5; + TIter Next(fPairs); + TObject *o = NULL; + while ((o=Next())) + fMap[col++] = fMatrix->AddColumn(TString(o->GetName())+".fDist"); +} + +void MEnergyEstParamDanielMkn421::StopMapping() +{ + fMatrix = NULL; + fPairs->Clear(); + fHillasSrc->Clear(); + fEnergy->Clear(); +} + +// -------------------------------------------------------------------------- +// +// Add a pair of input/output containers. +// eg. Add("MHillasSrc", "MEnergyEst"); +// Useful if you want to estimate the stuff for the source and antisource +// +void MEnergyEstParamDanielMkn421::Add(const TString hillas, const TString energy) +{ + fPairs->Add(new TNamed(hillas, energy)); + + AddToBranchList(hillas+".fDist"); +} + +void MEnergyEstParamDanielMkn421::Print(Option_t *opt) const +{ + for (int i=0; iSetEnergy(er); + est->SetImpact(ir); + est->SetReadyToSave(); + + } + + + return kTRUE; +} + + + + + + + + + + + + + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MEnergyEstParamDanielMkn421.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParamDanielMkn421.h (revision 2272) +++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParamDanielMkn421.h (revision 2272) @@ -0,0 +1,79 @@ +#ifndef MARS_MEnergyEstParamDanielMkn421 +#define MARS_MEnergyEstParamDanielMkn421 + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MHillasSrc; +class MEnergyEst; +class MHillas; +class MNewImagePar; +class MMcEvt; +class MHMatrix; + +class MEnergyEstParamDanielMkn421 : public MTask +{ +private: + Int_t fMap[100]; // FIXME! + + MHMatrix *fMatrix; + + TList *fPairs; + + TList *fHillasSrc; + TList *fEnergy; + + MHillas *fHillas; + MNewImagePar *fNewImagePar; + + TString fHillasName; + + TArrayD fA; + TArrayD fB; + + MMcEvt *fMc; + + void InitCoefficients(); + + Double_t GetVal(Int_t i) const; + Double_t fMm2Deg; + + Int_t PreProcess(MParList *plist); + Int_t Process(); + +public: + MEnergyEstParamDanielMkn421(const char *hil="MHillas", const char *name=NULL, const char *title=NULL); + ~MEnergyEstParamDanielMkn421(); + + void Add(const TString hillas, const TString energy="MEnergyEst"); + + void InitMapping(MHMatrix *mat); + void StopMapping(); + + Int_t GetNumCoeff() const { return fA.GetSize()+fB.GetSize(); } + + void SetCoeff(const TArrayD &arr); + void SetCoeffA(const TArrayD &arr); + void SetCoeffB(const TArrayD &arr); + + Double_t GetCoeff(Int_t i) { return i