Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6281) +++ trunk/MagicSoft/Mars/Changelog (revision 6282) @@ -52,4 +52,7 @@ common code from derived classes into the base-class - implemented CheckEnvLocal to handle the resource 'DataType' + - changed publicity of some variables and functions + to allow setting file type from outside + - added new function to set data-type for use in callisto.cc * mjobs/MJCalibTest.[h,cc], mjobs/MJCalibrateSignal.[h,cc], @@ -63,4 +66,13 @@ be able to write a function description about it (simple copy&past) + + * callisto.cc: + - added option for file-type + + * mjobs/MJOptimize.[h,cc], mjobs/MSequences.[h,cc], + mjobs/MJCut.[h,cc], ganymed.[cc,rc]: + - added to repository, but not yet to Makefile because + there is still some work to be done. But whoever is + interested in the new classes/program may already use it. Index: trunk/MagicSoft/Mars/callisto.cc =================================================================== --- trunk/MagicSoft/Mars/callisto.cc (revision 6281) +++ trunk/MagicSoft/Mars/callisto.cc (revision 6282) @@ -52,4 +52,8 @@ gLog << " -c Calculate the calibration constants" << endl; gLog << " -y Extract and calibrate signal" << endl << endl; + gLog << " Data Type (exclusive):" << endl; + gLog << " -raw Read input from raw-data" << endl; + gLog << " -mc Input root-files are monte carlo files" << endl; + gLog << " -root Read input from root-files (merpped) " << endl << endl; gLog << " Options:" << endl; gLog.Usage(); @@ -72,5 +76,6 @@ gLog << " --print-files Print Files taken from Sequence" << endl; gLog << " --print-only Do not execute anything except print" << endl; - gLog << " --use-test Apply cal. const. to the same cal. file (for testing)" << endl; + gLog << " --use-test Apply calibration constants to same calibration" << endl; + gLog << " file (for testing, calibration mode only)" << endl; gLog << " --config=callisto.rc Resource file [default=callisto.rc]" << endl; gLog << endl; @@ -141,4 +146,10 @@ Bool_t kModeY = arg.HasOnlyAndRemove("-y"); + MJCalib::DataType_t kDataType = MJCalib::kIsUseRootData; // root + if (arg.HasOnlyAndRemove("-raw")) + kDataType = MJCalib::kIsUseRawData; // raw + if (arg.HasOnlyAndRemove("-mc")) + kDataType = MJCalib::kIsUseMC; // monte carlo + if (!kInpathY.IsNull() || !kOutpathY.IsNull() || !kOutpath.IsNull() || !kPath.IsNull()) kModeY = kTRUE; @@ -146,12 +157,12 @@ kModeC = kTRUE; - if (!kModeC && !kModeY) - { - gLog << err << "Neither calibration (-c) nor signal extraction (-y) mode specified!" << endl; + if (!kModeC && !kModeY /*&& !kUseTest*/) + { + gLog << err << "Neither calibration (-c) nor signal extraction (-y) or test-mode (--use-test) specified!" << endl; Usage(); return 0; } - if (kModeC && kOutpathC.IsNull()) + if ((kModeC/* || kUseTest*/) && kOutpathC.IsNull()) kOutpathC = "."; if (kModeY) @@ -182,4 +193,5 @@ arg.Print("options"); gLog << endl; + return -1; } @@ -289,5 +301,5 @@ d->SetTitle(kSequence); - if (kModeC) + if (kModeC/* || kUseTest*/) { // @@ -302,4 +314,5 @@ job1.SetOverwrite(kOverwrite); job1.SetPathData(kInpathD); + job1.SetDataType(kDataType); job1.SetExtractionFundamental(); @@ -330,4 +343,5 @@ job2.SetOverwrite(kOverwrite); job2.SetPathData(kInpathD); + job2.SetDataType(kDataType); // job1.SetPathOut(kOutpathC); // not yet needed // job1.SetPathIn(kInpathC); // not yet needed @@ -351,31 +365,35 @@ } - // - // Do calibration - // - MJCalibration job3(Form("MJCalibration #%d", seq.GetSequence())); - job3.SetSequence(seq); - job3.SetEnv(kConfig); - job3.SetEnvDebug(kDebugEnv); - job3.SetDisplay(d);; - job3.SetOverwrite(kOverwrite); - job3.SetPathOut(kOutpathC); - job3.SetPathData(kInpathD); - // job2.SetPathIn(kInpathC); // not yet needed - - job3.SetBadPixels(job2.GetBadPixels()); - job3.SetExtractor(job2.GetExtractor()); - job3.SetExtractorCam(job2.GetPedestalCam()); - - if (!job3.ProcessFile(job1.GetPedestalCam())) - { - gLog << err << "Calculation of calibration failed." << endl << endl; - return -1; - } - - if (!job3.GetDisplay()) - { - gLog << warn << "Display closed by user... execution aborted." << endl << endl; - return 1; + if (kModeC) + { + // + // Do calibration + // + MJCalibration job3(Form("MJCalibration #%d", seq.GetSequence())); + job3.SetSequence(seq); + job3.SetEnv(kConfig); + job3.SetEnvDebug(kDebugEnv); + job3.SetDisplay(d); + job3.SetOverwrite(kOverwrite); + job3.SetPathOut(kOutpathC); + job3.SetPathData(kInpathD); + job3.SetDataType(kDataType); + // job2.SetPathIn(kInpathC); // not yet needed + + job3.SetBadPixels(job2.GetBadPixels()); + job3.SetExtractor(job2.GetExtractor()); + job3.SetExtractorCam(job2.GetPedestalCam()); + + if (!job3.ProcessFile(job1.GetPedestalCam())) + { + gLog << err << "Calculation of calibration failed." << endl << endl; + return -1; + } + + if (!job3.GetDisplay()) + { + gLog << warn << "Display closed by user... execution aborted." << endl << endl; + return 1; + } } @@ -383,5 +401,5 @@ { MJCalibTest job4(Form("MJCalibTest #%d", seq.GetSequence())); - job4.SetBadPixels(job3.GetBadPixels()); + job4.SetBadPixels(job2.GetBadPixels()); job4.SetSequence(seq); job4.SetEnv(kConfig); @@ -391,4 +409,5 @@ job4.SetPathOut(kOutpathC); job4.SetPathData(kInpathD); + job4.SetDataType(kDataType); if (!job4.ProcessFile(job1.GetPedestalCam())) @@ -423,4 +442,5 @@ job1.SetPathData(kInpathD); job1.SetPathIn(kInpathY); + job1.SetDataType(kDataType); //job1.SetPathOut(kOutpathY); // not yet needed job1.SetUseData(); @@ -452,4 +472,5 @@ job2.SetPathData(kInpathD); job2.SetPathIn(kInpathY); + job2.SetDataType(kDataType); // job1.SetPathOut(kOutpathC); // not yet needed // job1.SetPathIn(kInpathC); // not yet needed @@ -486,4 +507,5 @@ job3.SetPathData(kInpathD); job3.SetPathIn(kInpathY); + job3.SetDataType(kDataType); // job1.SetPathOut(kOutpathC); // not yet needed // job1.SetPathIn(kInpathC); // not yet needed @@ -519,4 +541,5 @@ job4.SetPathOut(kOutpathY); job4.SetPathData(kInpathD); + job4.SetDataType(kDataType); // Where to search for calibration files Index: trunk/MagicSoft/Mars/ganymed.cc =================================================================== --- trunk/MagicSoft/Mars/ganymed.cc (revision 6282) +++ trunk/MagicSoft/Mars/ganymed.cc (revision 6282) @@ -0,0 +1,230 @@ +#include +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MArgs.h" +#include "MArray.h" +#include "MDirIter.h" + +#include "MStatusDisplay.h" + +#include "MSequences.h" +#include "MJCut.h" + +using namespace std; + +static void StartUpMessage() +{ + gLog << all << endl; + + // 1 2 3 4 5 + // 12345678901234567890123456789012345678901234567890 + gLog << "========================================================" << endl; + gLog << " Ganymed - MARS V" << MARSVER << endl; + gLog << " MARS -- Gammas Are Now Your Most Exciting Discovery" << endl; + gLog << " Compiled on <" << __DATE__ << ">" << endl; + gLog << " Using ROOT v" << ROOTVER << endl; + gLog << "========================================================" << endl; + gLog << endl; +} + +static void Usage() +{ + // 1 2 3 4 5 6 7 8 + // 12345678901234567890123456789012345678901234567890123456789012345678901234567890 + gLog << all << endl; + gLog << "Sorry the usage is:" << endl; + gLog << " ganymed [-c] [-y] [options] sequences.txt" << endl << endl; + gLog << " Arguments:" << endl; + gLog << " sequences.txt: An ascii file defining a collection of sequences" << endl; + gLog << " Root Options:" << endl; + gLog << " -b Batch mode (no graphical output to screen)" << endl<0) + { + gLog << warn << "WARNING - Unknown commandline options..." << endl; + arg.Print("options"); + gLog << endl; + } + + // + // check for the right usage of the program + // + if (arg.GetNumArguments()!=1) + { + Usage(); + return -1; + } + + // + // Setup sequence file and check for its existance + // + const TString kSequences = arg.GetArgumentStr(0); + + if (gSystem->AccessPathName(kSequences, kFileExists)) + { + gLog << err << "Sorry, sequences file '" << kSequences << "' doesn't exist." << endl; + return -1; + } + + if (kDebugMem) + TObject::SetObjectStat(kTRUE); + + // + // Setup sequence and check its validity + // + MSequences seq(kSequences); + if (kNumAnalysis>=0) + seq.SetNumAnalysis(kNumAnalysis); + if (kPrintSeq || kPrintFiles) + { + gLog << all; + gLog.Separator(kSequences); + seq.Print(kPrintFiles?"files":""); + gLog << endl; + } + if (!seq.IsValid()) + { + gLog << err << "Sequences read but not valid!" << endl << endl; + return -1; + } + + // + // Initialize root + // + MArray::Class()->IgnoreTObjectStreamer(); + MParContainer::Class()->IgnoreTObjectStreamer(); + + TApplication app("Ganymed", &argc, argv); + if (!gROOT->IsBatch() && !gClient || gROOT->IsBatch() && !kBatch) + { + gLog << err << "Bombing... maybe your DISPLAY variable is not set correctly!" << endl; + return 1; + } + + // + // Update frequency by default = 1Hz + // + MStatusDisplay *d = new MStatusDisplay; + + // From now on each 'Exit' means: Terminate the application + d->SetBit(MStatusDisplay::kExitLoopOnExit); + d->SetTitle(kSequences); + + // + // Calculate pedestal for pedestal-calculation and calibration + // + MJCut job(Form("MJCut #%d", seq.GetNumAnalysis())); + job.SetEnv(kConfig); + job.SetEnvDebug(kDebugEnv); + job.SetDisplay(d);; + job.SetOverwrite(kOverwrite); + job.SetPathOut(kOutpath); + job.SetNameOutFile(kOutfile); + job.SetNameSummaryFile(kNameSummary); + job.SetNameResultFile(kNameResult); + job.EnableWriteOnly(kWriteOnly); + if (kWriteSummary) // Don't change flag set in SetNameSummaryFile + job.EnableStorageOfSummary(); + if (kWriteResult) // Don't change flag set in SetNameSummaryFile + job.EnableStorageOfResult(); + + if (!job.ProcessFile(seq)) + { + gLog << err << "Calculation of cuts failed." << endl << endl; + return -1; + } + + if (!job.GetDisplay()) + { + gLog << warn << "Display closed by user... execution aborted." << endl << endl; + return 1; + } + + if (kBatch || kQuit) + delete d; + else + { + // From now on each 'Close' means: Terminate the application + d->SetBit(MStatusDisplay::kExitLoopOnClose); + + // Wait until the user decides to exit the application + app.Run(kFALSE); + } + + if (TObject::GetObjectStat()) + { + TObject::SetObjectStat(kFALSE); + gObjectTable->Print(); + } + + return 0; +} Index: trunk/MagicSoft/Mars/ganymed.rc =================================================================== --- trunk/MagicSoft/Mars/ganymed.rc (revision 6282) +++ trunk/MagicSoft/Mars/ganymed.rc (revision 6282) @@ -0,0 +1,118 @@ +# ========================================================================== +############################################################################# +# ========================================================================== +# General +# ========================================================================== +############################################################################# +# ========================================================================== + +# ------------------------------------------------------------------------- +# Use this if you want to setup the logging stream for the jobs +# (overwrites command line options) +# ------------------------------------------------------------------------- +#MLog.VerbosityLevel: 2 +#MLog.DebugLevel: 1 +#MLog.NoColors: yes + +# ========================================================================== +############################################################################# +# ========================================================================== +# Cuts +# ========================================================================== +############################################################################# +# ========================================================================== + +# ------------------------------------------------------------------------- +# Use this if you want to write output to somewhere here you can give +# the output path +# ------------------------------------------------------------------------- +#MJCuts.PathOut: . + +# ------------------------------------------------------------------------- +# Use this to setup binnings. For more details see: MBinning::ReadEnv +# ------------------------------------------------------------------------- +#BinningAlpha.Raw: 9 0 90 +#BinningFalseSource.Raw: 30 -1.2 1.2 +#BinningEnergyEst.Raw: 25 10 1000000 log +#BinningTheta.Raw: 50 0 60 cos + +# ------------------------------------------------------------------------- +# Setup fit mode and ranges +# ------------------------------------------------------------------------- +#MAlphaFitter.SignalIntegralMax: 15 +#MAlphaFitter.SignalFitMax: 25 +#MAlphaFitter.BackgroundFitMin: 35 +#MAlphaFitter.BackgroundFitMax: 75 +#MAlphaFitter.ScaleMin: 35 +#MAlphaFitter.ScaleMax: 80 +#MAlphaFitter.PolynomOrder: 2 +#/*MISSING*/ ScaleMode + +# ------------------------------------------------------------------------- +# There are three cuts: +# Cut0: After energy estimation, before writing summary file +# Cut1: After writing summary file before filling false source plots +# Cut2: After filling false source plots before filling alpha plots +# Cut3: After filling alpha plots (eg. Alpha cuts) for image parameter +# display +# ------------------------------------------------------------------------- +Cut0.Inverted: Yes +Cut1.Inverted: Yes +Cut2.Inverted: Yes +Cut3.Inverted: Yes + +#Cut0.Condition: MCameraCooling.fTempCenter < 55 +Cut0.Condition: {0} && {1} +Cut0.0: log10(MNewImagePar.fConc) < 1.19 - (0.37 * log10(MHillas.fSize)) +Cut0.1: MNewImagePar.fConc > 0.45 + +# If you could setup MFEventSelector by resource file you could use it here +# To produce trainings and test sample use: "{MMcEvt.fEvtNumber%2}>0.5" +#Cut1.Condition: +#Cut1.Param0: 0.056 +#Cut1.Param8: 0.222 +#Cut1.Param16: 0.042 +#Cut1.Param24: 0.087 +#Cut1.Param32: 0.33 + +Cut1.Condition: {0} && {1} && {2} && {3} && {4} + +Cut1.0: MHillas.fLength*MGeomCam.fConvMm2Deg > 0.0063*log10(MHillas.fSize)+0.0974 +Cut1.1: MHillas.fLength*MGeomCam.fConvMm2Deg < 0.0040*log10(MHillas.fSize)+0.1791 +Cut1.2: MHillas.fWidth *MGeomCam.fConvMm2Deg > 0.0067*log10(MHillas.fSize)+0.0511 +Cut1.3: MHillas.fWidth *MGeomCam.fConvMm2Deg < 0.0050*log10(MHillas.fSize)+0.1203 +Cut1.4: MHillas.fSize>95 + +Cut2.Condition: MHillasSrc.fDist*MGeomCam.fConvMm2Deg>0.55 + +Cut3.Condition: abs(MHillasSrc.fAlpha)<10 + +# ------------------------------------------------------------------------- +# Setup cuts in false source plot +# ------------------------------------------------------------------------- +MHFalseSource.DistMin: 0.55 +#MHFalseSource.DistMax: 0.55 +#MHFalseSource.DWMin: 0.55 +#MHFalseSource.DWMax: 0.55 + +# ------------------------------------------------------------------------- +# Energy Estimation +# ------------------------------------------------------------------------- + +#EstimateEnergy: MEnergyEstimate +EstimateEnergy.Rule: {0} + {1} + {2} + {3} + {4} + {5} +EstimateEnergy.0: 6.3 +EstimateEnergy.1: 0.04*MHillasSrc.fDist*MGeomCam.fConvMm2Deg +EstimateEnergy.2: -0.13*MHillas.fLength*MGeomCam.fConvMm2Deg +EstimateEnergy.3: 0.15*MHillas.fSize +EstimateEnergy.4: 0.0000519*MHillas.fSize*MHillasSrc.fDist*MGeomCam.fConvMm2Deg +EstimateEnergy.5: 0.0000519*MHillas.fLength*MGeomCam.fConvMm2Deg*MHillasSrc.fDist*MGeomCam.fConvMm2Deg + + +# ------------------------------------------------------------------------- +# Use this to executa a task (eg to calc hadronness) after energy +# before Cut1 +# ------------------------------------------------------------------------- +# CalcHadronness: MRanForestCalc +# CalcHadronness.File: /home/tbretz/Mars.cvs/RFspot3cmAll.root +# CalcHadronness.Debug: No Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 6281) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 6282) @@ -140,4 +140,5 @@ // R use the range specified in the function range // Q quiet mode + // E Perform better Errors estimation using Minos technique h.Fit(fFunc, "NQI", "", bgmin, bgmax); @@ -175,4 +176,5 @@ // R use the range specified in the function range // Q quiet mode + // E Perform better Errors estimation using Minos technique h.Fit(fFunc, "NQI", "", 0, sigmax); @@ -261,4 +263,5 @@ // R use the range specified in the function range // Q quiet mode + // E Perform better Errors estimation using Minos technique TH1D h(hon); h.Add(&hof, -1); @@ -317,4 +320,5 @@ f.fScaleMode = fScaleMode; f.fScaleUser = fScaleUser; + f.fStrategy = fStrategy; f.fCoefficients.Set(fCoefficients.GetSize()); f.fCoefficients.Reset(); @@ -494,4 +498,16 @@ } +Double_t MAlphaFitter::GetMinimizationValue() const +{ + switch (fStrategy) + { + case kSignificance: + return -GetSignificance(); + case kSignificanceChi2: + return -GetSignificance()/GetChiSqSignal(); + } + return 0; +} + Int_t MAlphaFitter::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { @@ -537,4 +553,16 @@ } - return rc; -} + if (IsEnvDefined(env, prefix, "MinimizationStrategy", print)) + { + TString txt = GetEnvValue(env, prefix, "MinimizationStrategy", ""); + txt = txt.Strip(TString::kBoth); + txt.ToLower(); + if (txt==(TString)"Significance") + fStrategy = kSignificance; + if (txt==(TString)"SignificanceChi2") + fStrategy = kSignificanceChi2; + rc = kTRUE; + } + + return rc; +} Index: trunk/MagicSoft/Mars/mjobs/MJCalib.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCalib.h (revision 6281) +++ trunk/MagicSoft/Mars/mjobs/MJCalib.h (revision 6282) @@ -10,35 +10,36 @@ class MJCalib : public MJob { +public: + enum DataType_t + { + kIsUseRootData, + kIsUseRawData, + kIsUseMC + }; + private: - Byte_t fDataFlag; // Bit-field to store the data type + Byte_t fDataFlag; // Bit-field to store the data type protected: + MRunIter *fRuns; // Data files, only used for test applications, default is sequence files! - MRunIter *fRuns; // Data files, only used for test applications, default is sequence files! - - enum DataType_t - { - kIsUseRawData, - kIsUseRootData, - kIsUseMC - }; - - Bool_t IsUseRawData() const { return fDataFlag==kIsUseRawData; } - Bool_t IsUseRootData() const { return fDataFlag==kIsUseRootData; } - Bool_t IsUseMC() const { return fDataFlag==kIsUseMC; } - - void SetUseRawData ( const Bool_t b=kTRUE ) { fDataFlag=kIsUseRawData; } - void SetUseRootData( const Bool_t b=kTRUE ) { fDataFlag=kIsUseRootData; } - void SetUseMC ( const Bool_t b=kTRUE ) { fDataFlag=kIsUseMC; } - - Bool_t CheckEnvLocal(); + Bool_t CheckEnvLocal(); public: + MJCalib(); - MJCalib(); + void SetInput(MRunIter *iter) { fRuns = iter; } + + Bool_t IsUseRawData() const { return fDataFlag==kIsUseRawData; } + Bool_t IsUseRootData() const { return fDataFlag==kIsUseRootData; } + Bool_t IsUseMC() const { return fDataFlag==kIsUseMC; } + + void SetUseRawData () { fDataFlag=kIsUseRawData; } + void SetUseRootData() { fDataFlag=kIsUseRootData; } + void SetUseMC () { fDataFlag=kIsUseMC; } + + void SetDataType(DataType_t type) { fDataFlag=type; } - void SetInput ( MRunIter *iter ) { fRuns = iter; } - - ClassDef(MJCalib, 0) // Base class for calibration jobs + ClassDef(MJCalib, 0) // Base class for calibration jobs }; Index: trunk/MagicSoft/Mars/mjobs/MJCut.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 6282) +++ trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 6282) @@ -0,0 +1,490 @@ +/* ======================================================================== *\ +! +! * +! * 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, 1/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MJCut +// +// FIXME: Preparation for wobble mode missing +// +///////////////////////////////////////////////////////////////////////////// +#include "MJCut.h" + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MTaskList.h" +#include "MEvtLoop.h" + +#include "MStatusDisplay.h" + +#include "MReadReports.h" +#include "MPrint.h" +#include "MContinue.h" +#include "MEnergyEstimate.h" +#include "MTaskEnv.h" +#include "MSrcPosCalc.h" +#include "MHillasCalc.h" +#include "MFillH.h" +#include "MWriteRootFile.h" + +#include "../mhflux/MAlphaFitter.h" +#include "MBinning.h" +#include "MSequences.h" +#include "MParameters.h" +#include "MObservatory.h" + +ClassImp(MJCut); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MJCut::MJCut(const char *name, const char *title) + : fStoreSummary(kFALSE), fStoreResult(kFALSE), fWriteOnly(kFALSE), + fEstimateEnergy(0), fCalcHadronness(0) +{ + fName = name ? name : "MJCut"; + fTitle = title ? title : "Standard program to perform g/h-seperation cuts"; +} + +MJCut::~MJCut() +{ + if (fEstimateEnergy) + delete fEstimateEnergy; + if (fCalcHadronness) + delete fCalcHadronness; +} + +// -------------------------------------------------------------------------- +// +// Set the name of the summary file (events after cut0) +// If you give a name the storage of this file is enabled implicitly. +// If you give no filename the storage is neither enabled nor disabled, +// but the storage file name is reset. +// If no filename is set the default filename is used. +// You can explicitly enable or disable the storage using EnableStoreOf*() +// The default argument is no filename. +// +void MJCut::SetNameSummaryFile(const char *name) +{ + fNameSummary=name; + if (!fNameSummary.IsNull()) + EnableStorageOfSummary(); +} + +// -------------------------------------------------------------------------- +// +// Set the name of the summary file (events after cut3) +// If you give a name the storage of this file is enabled implicitly. +// If you give no filename the storage is neither enabled nor disabled, +// but the storage file name is reset. +// If no filename is set the default filename is used. +// You can explicitly enable or disable the storage using EnableStoreOf*() +// The default argument is no filename. +// +void MJCut::SetNameResultFile(const char *name) +{ + fNameResult=name; + if (!fNameResult.IsNull()) + EnableStorageOfResult(); +} + +// -------------------------------------------------------------------------- +// +// Setup a task estimating the energy. The given task is cloned. +// +void MJCut::SetEnergyEstimator(const MTask *task) +{ + if (fEstimateEnergy) + delete fEstimateEnergy; + fEstimateEnergy = task ? (MTask*)task->Clone() : 0; +} + +// -------------------------------------------------------------------------- +// +// Setup a task calculating the hadronness. The given task is cloned. +// +void MJCut::SetHadronnessCalculator(const MTask *task) +{ + if (fCalcHadronness) + delete fCalcHadronness; + fCalcHadronness = task ? (MTask*)task->Clone() : 0; +} + +// -------------------------------------------------------------------------- +// +// return fOutputPath+"/ganymed%08d.root", num +// +TString MJCut::GetOutputFile(UInt_t num) const +{ + TString p(fPathOut); + p += "/"; + p += fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data(); + return p; +} + +/* +Bool_t MJCut::ReadTasks(const char *fname, MTask* &env1, MTask* &env2) const +{ + // const TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence()); + + *fLog << inf << "Reading from file: " << fname << endl; + + TFile file(fname, "READ"); + if (!file.IsOpen()) + { + *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl; + return kFALSE; + } + + TObject *o = file.Get("EstimateEnergy"); + if (o && !o->InheritsFrom(MTask::Class())) + { + *fLog << err << dbginf << "ERROR - EstimateEnergy read from " << fname << " doesn't inherit from MTask!" << endl; + return kFALSE; + } + env1 = o ? (MTask*)o->Clone() : NULL; + + o = file.Get("CalcHadronness"); + if (o && !o->InheritsFrom(MTask::Class())) + { + *fLog << err << dbginf << "ERROR - CalcHadronness read from " << fname << " doesn't inherit from MTask!" << endl; + return kFALSE; + } + env2 = o ? (MTask*)o->Clone() : NULL; + + return kTRUE; +} +*/ + +// -------------------------------------------------------------------------- +// +// Write the tasks in cont to the file corresponding to analysis number num, +// see GetOutputFile() +// +Bool_t MJCut::WriteTasks(UInt_t num, TObjArray &cont) const +{ + if (fPathOut.IsNull()) + { + *fLog << inf << "No output path specified via SetPathOut - no output written." << endl; + return kTRUE; + } + + const TString oname(GetOutputFile(num)); + + *fLog << inf << "Writing to file: " << oname << endl; + + TFile file(oname, fOverwrite?"RECREATE":"NEW", "File created by MJCut", 9); + if (!file.IsOpen()) + { + *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl; + return kFALSE; + } + + return WriteContainer(cont); +} + +// -------------------------------------------------------------------------- +// +// Write the result plots and other results to the file corresponding to +// analysis number num, see GetOutputFile() +// +Bool_t MJCut::WriteResult(UInt_t num) const +{ + TObjArray arr; + return WriteContainer(arr, GetOutputFile(num), "UPDATE"); +} + +// -------------------------------------------------------------------------- +// +// MJCut allows to setup several option by a resource file: +// MJCut.WriteSummary: yes, no +// MJCut.SummaryFile: filename +// MJCut.WriteResult: yes, no +// MJCut.ResultFile: filename +// +Bool_t MJCut::CheckEnvLocal() +{ + const TString f0(GetEnv("SummaryFile", "")); + const TString f1(GetEnv("ResultFile", "")); + if (!f0.IsNull()) + SetNameSummaryFile(f0); + if (!f1.IsNull()) + SetNameResultFile(f1); + + EnableStorageOfSummary(GetEnv("SummaryFile", fStoreSummary)); + EnableStorageOfResult(GetEnv("ResultFile", fStoreResult)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Setup write to write: +// container tree optional? +// -------------- ---------- ----------- +// "MHillas" to "Events" +// "MHillasSrc" to "Events" +// "MHadronness" to "Events" yes +// "DataType" to "Events" +// +void MJCut::SetupWriter(MWriteRootFile &write) const +{ + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MHadronness", "Events", kFALSE); + write.AddContainer("DataType", "Events"); + + // Should not be the default: Either as option, or as + // setup from resource file + // write.AddContainer("MHillasExt", "Events"); + // write.AddContainer("MImagePar", "Events"); + // write.AddContainer("MNewImagePar", "Events"); +} + +Bool_t MJCut::ProcessFile(const MSequences &seq) +{ + if (!seq.IsValid()) + { + *fLog << err << "ERROR - Sequences invalid!" << endl; + return kFALSE; + } + + CheckEnv(); + + // -------------------------------------------------------------------------------- + + *fLog << inf; + fLog->Separator(GetDescriptor()); + *fLog << "Perform cuts for sequences " << seq.GetName() << endl; + *fLog << endl; + + // -------------------------------------------------------------------------------- + + // Setup Parlist + MParList plist; + plist.AddToList(this); // take care of fDisplay! + + MParameterI par("DataType"); + plist.AddToList(&par); + + // Setup Tasklist + MTaskList tlist; + plist.AddToList(&tlist); + + // La Palma Magic1 + MObservatory obs; + plist.AddToList(&obs); + + // Initialize binnings + MBinning bins1(18, 0, 90, "BinningAlpha", "lin"); + MBinning bins2(25, 10, 1000000, "BinningEnergyEst", "log"); + MBinning bins3(50, 0, 60, "BinningTheta", "cos"); + MBinning bins4("BinningFalseSource"); + plist.AddToList(&bins1); + plist.AddToList(&bins2); + plist.AddToList(&bins3); + plist.AddToList(&bins4); + + // -------------------------------------------------------------------------------- + + // Setup fitter and histograms + MAlphaFitter fit; + plist.AddToList(&fit); + + MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha"); + MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS"); + + // FIXME: If fPathIn read cuts and energy estimator from file! + MContinue cont0("", "Cut0"); + MContinue cont1("", "Cut1"); + MContinue cont2("", "Cut2"); + MContinue cont3("", "Cut3"); + cont0.SetAllowEmpty(); + cont1.SetAllowEmpty(); + cont2.SetAllowEmpty(); + cont3.SetAllowEmpty(); + + // ------------- Loop Off Data -------------------- + MReadReports readoff; + readoff.AddTree("Events", "MTime.", kTRUE); + readoff.AddTree("Drive"); + readoff.AddTree("EffectiveOnTime"); + seq.AddFilesOff(readoff); + + const TString path(Form("%s/", fPathOut.Data())); + TString fname0(path); + TString fname1(path); + fname0 += fNameSummary.IsNull() ? Form("ganymed%08d-summary.root", seq.GetNumAnalysis()) : fNameSummary; + fname1 += fNameResult.IsNull() ? Form("ganymed%08d-result.root", seq.GetNumAnalysis()) : fNameResult; + + MWriteRootFile write0(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW"); + MWriteRootFile write1(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW"); + if (CanStoreSummary()) + SetupWriter(write0); + if (CanStoreSummary()) + SetupWriter(write1); + + MEnergyEstimate est; + + MTaskEnv taskenv1("EstimateEnergy"); + taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est); + + MTaskEnv taskenv2("CalcHadronness"); + taskenv2.SetDefault(fCalcHadronness); + + MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre"); + MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost"); + fill1a.SetNameTab("PreCut"); + fill2a.SetNameTab("PostCut"); + + MPrint print2("MEffectiveOnTime"); + + // How to get source position from off- and on-data? + MSrcPosCalc scalc; + + //MHillasCalc hcalc; + //hcalc.SetFlags(MHillasCalc::kCalcHillasSrc); + + MTaskList tlist2; + tlist2.AddToList(&taskenv1); + tlist2.AddToList(&taskenv2); + tlist2.AddToList(&cont0); + if (CanStoreSummary()) + tlist2.AddToList(&write0); + if (!fWriteOnly) + tlist2.AddToList(&fill1a); + tlist2.AddToList(&cont1); + if (!fWriteOnly) + tlist2.AddToList(&ffs); + //tlist2.AddToList(&scalc); + //tlist2.AddToList(&hcalc); + tlist2.AddToList(&cont2); + if (!fWriteOnly) + tlist2.AddToList(&fill2a); + if (!fWriteOnly) + tlist2.AddToList(&falpha); + tlist2.AddToList(&cont3); + if (CanStoreResult()) + tlist2.AddToList(&write1); + + tlist.AddToList(&readoff); + tlist.AddToList(&print2, "EffectiveOnTime"); + tlist.AddToList(&tlist2, "Events"); + + par.SetVal(0); + + // Create and setup the eventloop + MEvtLoop evtloop(fName); + evtloop.SetParList(&plist); + evtloop.SetDisplay(fDisplay); + evtloop.SetLogStream(fLog); + if (!SetupEnv(evtloop)) + return kFALSE; + + if (seq.HasOffSequences()) + { + // Execute first analysis + if (!evtloop.Eventloop(fMaxEvents)) + { + *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl; + return kFALSE; + } + + tlist.PrintStatistics(); + + if (!evtloop.GetDisplay()) + { + *fLog << err << GetDescriptor() << ": Execution stopped by used." << endl; + return kFALSE; + } + } + + // ------------- Loop On Data -------------------- + MReadReports readon; + readon.AddTree("Events", "MTime.", kTRUE); + readon.AddTree("Drive"); + readon.AddTree("EffectiveOnTime"); + seq.AddFilesOn(readon); + + MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre"); + MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost"); + fill1b.SetNameTab("PreCut"); + fill2b.SetNameTab("PostCut"); + fill1b.SetDrawOption("same"); + fill2b.SetDrawOption("same"); + + MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha"); + MFillH ffs2("MHFalseSource", "MHillas", "FillFS"); + + tlist.Replace(&readon); + if (!fWriteOnly) + { + tlist2.Replace(&fill1b); + tlist2.Replace(&fill2b); + tlist2.Replace(&falpha2); + tlist2.Replace(&ffs2); + } + + par.SetVal(1); + + TObjArray cont; + cont.Add(&cont0); + cont.Add(&cont1); + cont.Add(&cont2); + cont.Add(&cont3); + if (taskenv1.GetTask()) + cont.Add(taskenv1.GetTask()); + if (taskenv2.GetTask()) + cont.Add(taskenv2.GetTask()); + + if (!WriteTasks(seq.GetNumAnalysis(), cont)) + return kFALSE; + + // Execute first analysis + if (!evtloop.Eventloop(fMaxEvents)) + { + *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl; + return kFALSE; + } + + tlist.PrintStatistics(); + + // FIXME: Perform fit and plot energy dependant alpha plots + // and fit result to new tabs! + if (!WriteResult(seq.GetNumAnalysis())) + return kFALSE; + + *fLog << all << GetDescriptor() << ": Done." << endl; + *fLog << endl << endl; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mjobs/MJCut.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCut.h (revision 6282) +++ trunk/MagicSoft/Mars/mjobs/MJCut.h (revision 6282) @@ -0,0 +1,55 @@ +#ifndef MARS_MJCut +#define MARS_MJCut + +#ifndef MARS_MJob +#include "MJob.h" +#endif + +class MTask; +class MSequences; +class MWriteRootFile; + +class MJCut : public MJob +{ +private: + Bool_t fStoreSummary; + Bool_t fStoreResult; + Bool_t fWriteOnly; + + TString fNameSummary; + TString fNameResult; + TString fNameOutput; + + MTask *fEstimateEnergy; + MTask *fCalcHadronness; + + TString GetOutputFile(UInt_t num) const; + Bool_t CheckEnvLocal(); + void SetupWriter(MWriteRootFile &write) const; + Bool_t WriteTasks(UInt_t num, TObjArray &cont) const; + Bool_t WriteResult(UInt_t num) const; + + Bool_t CanStoreSummary() const { return !fPathOut.IsNull() && fStoreSummary; } + Bool_t CanStoreResult() const { return !fPathOut.IsNull() && fStoreResult; } + +public: + MJCut(const char *name=NULL, const char *title=NULL); + ~MJCut(); + + Bool_t ProcessFile(const MSequences &seq); + + void EnableStorageOfSummary(Bool_t b=kTRUE) { fStoreSummary = b; } // See SetNameSummary + void EnableStorageOfResult(Bool_t b=kTRUE) { fStoreResult = b; } // See SetNameResult + void EnableWriteOnly(Bool_t b=kTRUE) { fWriteOnly = b; } + + void SetNameSummaryFile(const char *name=""); + void SetNameResultFile(const char *name=""); + void SetNameOutFile(const char *name="") { fNameOutput=name; } + + void SetEnergyEstimator(const MTask *task=0); + void SetHadronnessCalculator(const MTask *task=0); + + ClassDef(MJCut, 0) // Standard program to perform g/h-seperation cuts +}; + +#endif Index: trunk/MagicSoft/Mars/mjobs/MJOptimize.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimize.cc (revision 6282) +++ trunk/MagicSoft/Mars/mjobs/MJOptimize.cc (revision 6282) @@ -0,0 +1,1010 @@ +/* ======================================================================== *\ +! +! * +! * 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/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MJOptimize +// +// Class for otimizing the parameters of the supercuts +// +// Minimization Control +// ==================== +// +// To choose the minimization algorithm use: +// void SetOptimizer(Optimizer_t o); +// +// Allowed options are: +// enum Optimizer_t +// { +// kMigrad, // Minimize by the method of Migrad +// kSimplex, // Minimize by the method of Simplex +// kMinimize, // Migrad + Simplex (if Migrad fails) +// kMinos, // Minos error determination +// kImprove, // Local minimum search +// kSeek, // Minimize by the method of Monte Carlo +// kNone // Skip optimization +// }; +// +// For more details on the methods see TMinuit. +// +// +// You can change the behaviour of the minimization using +// +// void SetNumMaxCalls(UInt_t num=0); +// void SetTolerance(Float_t tol=0); +// +// While NumMaxCalls is the first, Tolerance the second arguement. +// For more details start root and type +// +// gMinuit->mnhelp("command") +// +// while command can be +// * MIGRAD +// * SIMPLEX +// * MINIMIZE +// * MINOS +// * IMPROVE +// * SEEK +// +// The default (num==0 and tol==0) should always give you the +// corresponding defaults used in Minuit. +// +// +// FIXME: Implement changing cut in hadronness... +// FIXME: Show MHSignificance on MStatusDisplay during filling... +// FIXME: Choose step-size percentage as static data membewr +// FIXME: Choose minimization method +// +///////////////////////////////////////////////////////////////////////////// +#include "MJOptimize.h" + +#include +#include + +#include + +#include + +#include +#include + +#include "MParList.h" +#include "MTaskList.h" +#include "MGeomCamCT1.h" +#include "MFEventSelector.h" +#include "MReadTree.h" +#include "MHMatrix.h" +#include "MEnergyEstimate.h" +#include "MMatrixLoop.h" +#include "MChisqEval.h" +#include "MEvtLoop.h" +#include "MDataElement.h" +#include "MDataMember.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MParameters.h" +#include "MFillH.h" +#include "MF.h" +#include "MFilterList.h" +#include "../mfilter/MFSupercuts.h" +#include "MContinue.h" +#include "MGeomCamMagic.h" +#include "../mimage/MHillasSrcCalc.h" +#include "MHMatrix.h" +#include "MMatrixLoop.h" +#include "../mhflux/MHAlpha.h" +#include "MStatusDisplay.h" +#include "../mhflux/MHEnergyEst.h" +#include "MDirIter.h" +#include "../mjobs/MSequence.h" + +using namespace std; + +//------------------------------------------------------------------------ +// +// fcn calculates the function to be minimized (using TMinuit::Migrad) +// +void MJOptimize::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + MJOptimize *optim = (MJOptimize*)gMinuit->GetObjectFit(); + + // WORKAROUND --- FOR WHAT? + if (gMinuit->fEpsi<1e-2) + { + *optim->fLog << warn << "WARNING - For unknown reasons: fEspi<1e-100... resetting to 0.01." << endl; + gMinuit->fEpsi = 0.01; + } + + TMinuit *minuit = gMinuit; + f = optim->Fcn(TArrayD(TMath::Min(gMinuit->fMaxpar, optim->fParameters.GetSize()), par), minuit); + gMinuit = minuit; + +} + +Double_t MJOptimize::Fcn(const TArrayD &par, TMinuit *minuit) +{ + if (fEvtLoop->GetDisplay()!=fDisplay) + return 0; + /* + switch(iflag) + { + case 1: // first call + case 2: // calc derivative + break; + case 3: + // last call + MStatusDisplay *d = new MStatusDisplay; + fEvtLoop->SetDisplay(d); + break; + } + */ + MParList *plist = fEvtLoop->GetParList(); + + MParameterD *eval = (MParameterD*)plist->FindObject("MinimizationValue", "MParameterD"); + MParContainer *pars = (MParContainer*)plist->FindObject("MParameters", "MParContainer"); + + MRead *read = (MRead*)plist->FindObject("MTaskList")->FindObject("MRead"); + if (read) + read->Rewind(); + + if (fDebug>=0) + { + *fLog << inf << "New Set: "; + for (Int_t i=0; iSetVariables(par); + eval->SetVal(0); + + if (fDebug<3) + gLog.SetNullOutput(kTRUE); + + TStopwatch clock; + clock.Start(); + fEvtLoop->Eventloop(fNumEvents); + clock.Stop(); + + if (fDebug<3) + gLog.SetNullOutput(kFALSE); + + const Double_t f = eval->GetVal(); + + if (fDebug>=0) + *fLog << inf << "Result F=" << f << endl; + + if (fDebug>=1 && minuit) + { + Double_t fmin, fedm, errdef; + Int_t n1, n2, istat; + minuit->mnstat(fmin, fedm, errdef, n1, n2, istat); + *fLog << inf << underline << "Minimization Status so far:" << endl; + *fLog << " Calls: " << minuit->fNfcn << " (max=" << gMinuit->fMaxIterations << ")" << endl; + *fLog << " Parameters: fixed=" << gMinuit->fNpfix << ", free=" << gMinuit->fNpar << endl; + *fLog << " Func min: " << fmin << " (Epsi=" << gMinuit->fEpsi << ", Apsi=" << gMinuit->fApsi << ")" << endl; + *fLog << " Found edm: " << fedm << endl; + *fLog << " ErrDef: " << errdef << endl; + *fLog << " Status: "; + + switch (istat) + { + case 0: *fLog << "n/a" << endl; break; + case 1: *fLog << "approximation only, not accurate" << endl; break; + case 2: *fLog << "full matrix, but forced positive-definite" << endl; break; + case 3: *fLog << "full accurate covariance matrix" << endl; break; + default: *fLog << "undefined" << endl; break; + } + } + + if (fDebug>=1) + { + clock.Print(); + fEvtLoop->GetTaskList()->PrintStatistics(); + } + + return f; +} + +MJOptimize::MJOptimize() : fDebug(-1), fNumEvents(0), fType(kSimplex), fNumMaxCalls(0), fTolerance(0) +{ + fRules.SetOwner(); + fFilter.SetOwner(); + + fNamesOn.SetOwner(); + fNamesOff.SetOwner(); +} + +//------------------------------------------------------------------------ +// +// Add seqeunces from list to reader +// +Bool_t MJOptimize::AddSequences(MRead &read, TList &list) const +{ + MDirIter files; + + TIter Next(&list); + TObject *o=0; + while ((o=Next())) + { + MSequence seq(o->GetName()); + if (!seq.IsValid()) + return kFALSE; + + seq.SetupDatRuns(files, o->GetTitle(), "I"); + } + + return read.AddFiles(files)>0; +} + +//------------------------------------------------------------------------ +// +// Add on-sequences: +// - fname: sequence file name (with path) +// - dir: directory were image files are stored +// +void MJOptimize::AddSequenceOn(const char *fname, const char *dir) +{ + fNamesOn.Add(new TNamed(fname, dir)); +} + +//------------------------------------------------------------------------ +// +// Add off-sequences: +// - fname: sequence file name (with path) +// - dir: directory were image files are stored +// +void MJOptimize::AddSequenceOff(const char *fname, const char *dir) +{ + fNamesOff.Add(new TNamed(fname, dir)); +} + +//------------------------------------------------------------------------ +// +// Empty list of on- and off-sequences +// +void MJOptimize::ResetSequences() +{ + fNamesOn.Delete(); + fNamesOff.Delete(); +} + +//------------------------------------------------------------------------ +// +// Add a parameter used in your filters (see AddFilter) The parameter +// index is returned, +// +// Int_t idx = AddParameter("log10(MHillas.fSize)"); +// +// The indices area starting with 0 always. +// +Int_t MJOptimize::AddParameter(const char *rule) +{ + fRules.Add(new TNamed(rule, "")); + return fRules.GetSize()-1; +} + +//------------------------------------------------------------------------ +// +// Add a filter which can be applied in the optimization (for deatils +// see correspodning Run function) You can use the indices you got by +// AddParameter, eg +// +// AddFilter("M[0] < 3.2"); +// +// if used in optimization you can do +// +// AddFilter("M[0] < [0]"); +// +// for more details, see SetParameter and FixParameter +// +void MJOptimize::AddFilter(const char *rule) +{ + fFilter.Add(new MF(rule)); +} + +//------------------------------------------------------------------------ +// +// Add a cut which is used to fill the matrix, eg "MMcEvt.fOartId<1.5" +// (The rule is applied, nit inverted: The matrix is filled with +// the events fullfilling the condition) +// +void MJOptimize::AddPreCut(const char *rule) +{ + fPreCuts.Add(new MF(rule)); +} + +//------------------------------------------------------------------------ +// +// Set the fParameters Array accoring to par. +// +void MJOptimize::SetParameters(const TArrayD &par) +{ + fParameters = par; +} + +//------------------------------------------------------------------------ +// +// Set the number of events processed by the eventloop. (Be carfull, +// if you are doing on-off analysis and you only process the first +// 1000 events which are on-events only the optimization may not work) +// +void MJOptimize::SetNumEvents(UInt_t n) +{ + fNumEvents=n; +} + +//------------------------------------------------------------------------ +// +// Set a debug level, which tells the optimization how much information +// is displayed about and in the running eventloop. +// +void MJOptimize::SetDebug(UInt_t n) +{ + fDebug=n; +} + +//------------------------------------------------------------------------ +// +// Set a optimization algorithm to be used. For more information see +// TMinuit. +// +// Available Algorithms are: +// kMigrad, // Minimize by the method of Migrad +// kSimplex, // Minimize by the method of Simplex +// kSeek // Minimize by the method of Monte Carlo +// +void MJOptimize::SetOptimizer(Optimizer_t o) +{ + fType = o; +} + +//------------------------------------------------------------------------ +// +// If a status didplay is set, search for tab "Optimizer". +// If not found, create it. +// In the tab search for TMultiGraph "Parameters". +// If not found create it. +// If empty create TGraphs. +// Check number of graphs vs. number of parameters. +// return TList with graphs. +// +TList *MJOptimize::GetPlots() const +{ + if (!fDisplay) + return NULL; + + TCanvas *c=fDisplay->GetCanvas("Optimizer"); + if (!c) + c = &fDisplay->AddTab("Optimizer"); + + TMultiGraph *mg = dynamic_cast(c->FindObject("Parameters")); + if (!mg) + mg = new TMultiGraph("Parameters", "Parameters of optimization"); + + TList *l = mg->GetListOfGraphs(); + if (!l) + { + const Int_t n = fParameters.GetSize(); + for (int i=0; iSetLineColor(kBlue); + mg->Add(g, ""); + AddPoint(mg->GetListOfGraphs(), i, i==n?1:fParameters[i]); + } + mg->SetBit(kCanDelete); + mg->Draw("al*"); + + l = mg->GetListOfGraphs(); + } + + return l->GetSize() == fParameters.GetSize()+1 ? l : NULL; +} + +//------------------------------------------------------------------------ +// +// Add a point with y=val as last point in idx-th Tgraph of list l. +// +void MJOptimize::AddPoint(TList *l, Int_t idx, Float_t val) const +{ + if (!l) + return; + + TGraph *gr = (TGraph*)l->At(idx); + gr->SetPoint(gr->GetN(), gr->GetN(), val); +} + +Int_t MJOptimize::Minuit(TMinuit &minuit, const char *cmd) const +{ + Int_t err; + Double_t tmp[2] = { fNumMaxCalls, fTolerance }; + minuit.mnexcm(cmd, tmp, 2, err); + + switch (err) + { + case 0: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm excuted normally." << endl; + break; + case 1: + *fLog << warn << GetDescriptor() << " TMinuit::mnexcm command is blank... ignored." << endl; + break; + case 2: + *fLog << warn << GetDescriptor() << " TMinuit::mnexcm command-line syntax error... ignored." << endl; + break; + case 3: + *fLog << warn << GetDescriptor() << " TMinuit::mnexcm unknown command... ignored." << endl; + break; + case 4: + *fLog << warn << GetDescriptor() << " TMinuit::mnexcm - Abnormal termination (eg Migrad not converged)" << endl; + break; + /* + case 5: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - Parameters requested." << endl; + break; + case 6: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET INPUT returned." << endl; + break; + case 7: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET TITLE returned." << endl; + break; + case 8: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET COVAR returned." << endl; + break; + case 9: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - reserved." << endl; + break; + case 10: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - END returned." << endl; + break; + case 11: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - EXIT or STOP returned." << endl; + break; + case 12: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - RETURN returned." << endl; + break;*/ + } + + return err; +} + +Bool_t MJOptimize::Optimize(MEvtLoop &evtloop) +{ + if (fType==kNone) + return kTRUE; + + gMinuit = new TMinuit(fParameters.GetSize()); + //gMinuit->SetPrintLevel(-1); + + gMinuit->SetFCN(fcn); + gMinuit->SetObjectFit(this); + + // + // Set starting values and step sizes for parameters + // + for (Int_t i=0; iDefineParameter(i, name, vinit, step, limlo, limup); + if (rc) + { + *fLog << err << dbginf << "Error in defining parameter #" << i << endl; + return kFALSE; + } + + if (step==0) + gMinuit->FixParameter(i); + } + + fEvtLoop = &evtloop; + + TList *g=GetPlots(); + + // Now ready for minimization step: + TStopwatch clock; + clock.Start(); + switch (fType) + { + case kSimplex: + Simplex(*gMinuit); + break; + case kMigrad: + Migrad(*gMinuit); + break; + case kMinimize: + Minimize(*gMinuit); + break; + case kMinos: + Minos(*gMinuit); + break; + case kImprove: + Improve(*gMinuit); + break; + case kSeek: + Seek(*gMinuit); + break; + case kNone: // Should never happen + return kFALSE; + } + clock.Stop(); + clock.Print(); + + if (evtloop.GetDisplay()!=fDisplay) + { + *fLog << inf << "Optimization aborted by user." << endl; + fDisplay = 0; + return kFALSE; + } + + *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl; + + // + // Update values of fA, fB: + // + for (Int_t i=0; iGetParameter(i,x1,x2); + fParameters[i] = x1; + cout << i << ": " << fParameters[i] << endl; + + AddPoint(g, i, x1); + } + AddPoint(g, fParameters.GetSize(), gMinuit->fAmin); + + delete gMinuit; + + return kTRUE; +} + +//------------------------------------------------------------------------ +// +// Optimize allows to use the optimizing by an eventloop based on +// some requirements. +// +// 1) The tasklist to be executed must have the name MTaskList and +// be an entry in the parameterlist. +// +// 2) The reading task (MReadMarsFile, MMatrixLoop) must have the name +// "MRead". If it derives from MRead Rewind() must be implemented, +// otherwise it must start reading from scratch each time its +// PreProcess is called. +// +// 3) The parameters to be optimized must be accessed through (currently) +// a single parameter container (MParContainer) called MParameters. +// The parameters are set through SetVariables. +// +// 4) The result of a single function call for minimization (eg. chisquare) +// must be available after the eventloop in a container of type +// MParameterD with the name "MinimizationResult". +// +// 5) The parameters to start with must have been set using +// MJOptimize::SetParameter or MJOptimize::SetParameters and +// MJOptimize::FixParameter +// +// The behaviour of the optimization can be changed using: +// void SetNumEvents(UInt_t n); +// void SetDebug(UInt_t n); +// void SetOptimizer(Optimizer_t o); +// +// After optimization the resulting parameters are set and another eventloop +// with a MStatusDisplay is set is called. The resulting parameters can be +// accessed using: GetParameters() +// +// To be fixed: +// - MStatusDisplay should show status while optimization is running +// - write result into MStatusDisplay +// - save result +// +Bool_t MJOptimize::Optimize(MParList &parlist) +{ + // Checks to make sure, that fcn doesn't crash + if (!parlist.FindCreateObj("MParameterD", "MinimizationValue")) + return kFALSE; + + if (!parlist.FindObject("MParameters", "MParContainer")) + { + *fLog << err << "MParameters [MParContainer] not found... abort." << endl; + return kFALSE; + } + + TString txt("Starting "); + switch (fType) + { + case kMigrad: txt += "Migrad"; break; + case kMinimize: txt += "Minimize"; break; + case kMinos: txt += "Minot"; break; + case kImprove: txt += "Improve"; break; + case kSimplex: txt += "Simplex"; break; + case kSeek: txt += "Seek"; break; + case kNone: txt += "no"; break; + } + txt += " optimization"; + + fLog->Separator(txt); + + // Setup eventloop + MEvtLoop evtloop; + evtloop.SetParList(&parlist); + evtloop.SetDisplay(fDisplay); // set display for evtloop and all childs + parlist.SetDisplay(0); // reset display for all contents of parlist and tasklist + evtloop.SetPrivateDisplay(); // prevent display from being cascaded again in PreProcess + + *fLog << inf << "Number of Parameters: " << fParameters.GetSize() << endl; + + if (!Optimize(evtloop)) + return kFALSE; + + gMinuit = 0; + + fEvtLoop->SetDisplay(fDisplay); + return Fcn(fParameters); + // Done already in Fcn + // list.SetVariables(fParameters); +} + +void MJOptimize::AddRulesToMatrix(MHMatrix &m) const +{ + TIter Next1(&fRules); + TObject *o1=0; + while ((o1=Next1())) + m.AddColumn(o1->GetName()); +} + +//------------------------------------------------------------------------ +// +// Fill matrix m by using read. Use rules as a filter if userules. +// +Bool_t MJOptimize::FillMatrix(MReadTree &read, MParList &parlist, Bool_t userules) +{ + MHMatrix *m = (MHMatrix*)parlist.FindObject("M", "MHMatrix"); + if (!m) + { + *fLog << err << "MJOptimize::FillMatrix - ERROR: M [MHMatrix] not found in parlist... abort." << endl; + return kFALSE; + } + + m->Print("cols"); + + //MParList parlist; + + // MGeomCamMagic cam; + // parlist.AddToList(&cam); + + MTaskList tlist; + parlist.Replace(&tlist); + + MFillH fillh(m); + + tlist.AddToList(&read); + + MFilterList list; + if (!list.AddToList(fPreCuts)) + *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl; + if (userules) + SetupFilters(list); + list.SetName("PreCuts"); // reset Name set by SetupFilters + list.SetInverted(kFALSE); // reset inversion set by SetupFilters + fillh.SetFilter(&list); + tlist.AddToList(&list); + + tlist.AddToList(&fillh); + + MEvtLoop fillloop; + fillloop.SetParList(&parlist); + fillloop.SetDisplay(fDisplay); + if (!fillloop.Eventloop(fNumEvents)) + { + *fLog << err << "Filling matrix failed..." << endl; + return kFALSE; + } + tlist.PrintStatistics(); + + *fLog << inf << "Read events from file '" << read.GetFileName() << "'" << endl; + + if (fillloop.GetDisplay()!=fDisplay) + { + fDisplay = 0; + *fLog << inf << "Optimization aborted by user." << endl; + return kFALSE; + } + + m->Print("size"); + + return kTRUE; +} + +//------------------------------------------------------------------------ +// +// Adds all filters to MFilterList +// +void MJOptimize::SetupFilters(MFilterList &list, MFilter *filter) const +{ + list.SetName("MParameters"); + list.SetInverted(); + + if (filter) + { + if (fFilter.GetSize()>0) + { + *fLog << inf; + *fLog << "INFORMATION - You are using an external filter and internal filters." << endl; + *fLog << " Please make sure that all parameters '[i]' are starting" << endl; + *fLog << " behind the number of parameters of the external filter." << endl; + } + list.AddToList(filter); + } + + if (!list.AddToList(fFilter)) + *fLog << err << "ERROR - Calling MFilterList::AddToList fFilter failed!" << endl; + + *fLog << inf << "Filter: "; + list.Print(); + *fLog << endl; +} + +//------------------------------------------------------------------------ +// +Bool_t MJOptimize::Run(const char *fname, MFilter *filter, MAlphaFitter *fit) +{ + if (fParameters.GetSize()==0) + { + *fLog << err << "Sorry, no parameters defined." << endl; + return kFALSE; + } + + fLog->Separator("Preparing On-only-optimization"); + + MParList parlist; + + MGeomCamMagic geom; // For GetConvMm2Deg + parlist.AddToList(&geom); + + MHMatrix m("M"); + AddRulesToMatrix(m); + parlist.AddToList(&m); + + MHAlpha hist; + hist.SkipHistTime(); + hist.SkipHistTheta(); + hist.SkipHistEnergy(); + hist.InitMapping(&m); + + if (filter && filter->InheritsFrom(MFSupercuts::Class())) + ((MFSupercuts*)filter)->InitMapping(&m); + + MReadTree read("Events"); + read.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! + if (fname) + read.AddFile(fname); + else + AddSequences(read, fNamesOn); + if (!FillMatrix(read, parlist)) + return kFALSE; + + MTaskList tasklist; + parlist.Replace(&tasklist); + if (fit) + parlist.AddToList(fit); + + MFilterList list; + SetupFilters(list, filter); + + MContinue contin(&list); + parlist.AddToList(&list); + + MFillH fill(&hist); + + MMatrixLoop loop(&m); + + tasklist.AddToList(&loop); + tasklist.AddToList(&list); + tasklist.AddToList(&contin); + tasklist.AddToList(&fill); + + if (!Optimize(parlist)) + return kFALSE; + + TObjArray cont; + cont.Add(&contin); + return WriteContainer(cont, fNameOut); +} + +//------------------------------------------------------------------------ +// +// Make sure, that filter->GetDataMember is correctly implemented!!!! +// +Bool_t MJOptimize::RunOnOff(const char *fname, MFilter *filter, MAlphaFitter *fit, const char *tree) +{ + if (fParameters.GetSize()==0) + { + *fLog << err << "Sorry, no parameters defined." << endl; + return kFALSE; + } + + fLog->Separator("Preparing On/Off-optimization"); + + MParList parlist; + + MGeomCamMagic geom; // For GetConvMm2Deg + parlist.AddToList(&geom); + + MHMatrix m("M"); + AddRulesToMatrix(m); + parlist.AddToList(&m); + + const Int_t idxdatatype = m.AddColumn("DataType.fVal"); + + MHAlpha histon, histof("MHAlphaOff"); + histon.SkipHistTime(); + histon.SkipHistTheta(); + //histon.SkipHistEnergy(); + histof.SkipHistTime(); + histof.SkipHistTheta(); + //histof.SkipHistEnergy(); + // FIXME: MHillasSrc.fAlpha is added twice! + histon.InitMapping(&m, 1); + histof.InitMapping(&m, 1); + + if (filter && filter->InheritsFrom(MFSupercuts::Class())) + ((MFSupercuts*)filter)->InitMapping(&m); + + parlist.AddToList(&histon); + parlist.AddToList(&histof); + + if (fname) + { + MReadTree read(tree); + read.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! + read.AddFile(fname); + if (!FillMatrix(read, parlist)) + return kFALSE; + } + else + { + MParameterI par("DataType"); + parlist.AddToList(&par); + + gLog.Separator("Reading On-Data"); + par.SetVal(1); + MReadTree readon(tree); + readon.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! + AddSequences(readon, fNamesOn); + if (!FillMatrix(readon, parlist)) + return kFALSE; + + gLog.Separator("Reading Off-Data"); + par.SetVal(0); + MReadTree readoff(tree); + readoff.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! + AddSequences(readoff, fNamesOff); + if (!FillMatrix(readoff, parlist)) + return kFALSE; + } + + MTaskList tasklist; + parlist.Replace(&tasklist); + if (fit) + parlist.AddToList(fit); + + MFilterList list; + SetupFilters(list, filter); + + MContinue contin(&list); + parlist.AddToList(&list); + + MFillH fillon(&histon, "", "FillHistOn"); + MFillH fillof(&histof, "", "FillHistOff"); + + MF f0(Form("M[%d]<0.5", idxdatatype), "FilterOffData"); + MF f1(Form("M[%d]>0.5", idxdatatype), "FilterOnData"); + + fillof.SetFilter(&f0); + fillon.SetFilter(&f1); + + MMatrixLoop loop(&m); + + tasklist.AddToList(&loop); + tasklist.AddToList(&list); + tasklist.AddToList(&contin); + tasklist.AddToList(&f0); + tasklist.AddToList(&f1); + tasklist.AddToList(&fillof); + tasklist.AddToList(&fillon); + + if (!Optimize(parlist)) + return kFALSE; + + TObjArray cont; + cont.Add(&contin); + return WriteContainer(cont, fNameOut); +} + +//------------------------------------------------------------------------ +// +// Read all events from file which do match rules and optimize +// energy estimator. +// +Bool_t MJOptimize::RunEnergy(const char *fname, const char *rule) +{ + if (fParameters.GetSize()==0) + { + *fLog << err << "Sorry, no parameters defined." << endl; + return kFALSE; + } + + fLog->Separator("Preparing Energy optimization"); + + MParList parlist; + + MGeomCamMagic geom; // For GetConvMm2Deg + parlist.AddToList(&geom); + + MHMatrix m("M"); + AddRulesToMatrix(m); + parlist.AddToList(&m); + + MEnergyEstimate est("MParameters"); + est.SetRule(rule); + parlist.AddToList(&est); + + MHEnergyEst hist; + hist.InitMapping(&m); + + MReadTree read("Events"); + //read.DisableAutoScheme(); + if (fname) + read.AddFile(fname); + else + AddSequences(read, fNamesOn); + if (!FillMatrix(read, parlist, kTRUE)) + return kFALSE; + + MTaskList tasklist; + parlist.Replace(&tasklist); + + MFillH fill(&hist); + + MMatrixLoop loop(&m); + + tasklist.AddToList(&loop); + tasklist.AddToList(&est); + tasklist.AddToList(&fill); + + if (!Optimize(parlist)) + return kFALSE; + + TObjArray cont; + cont.Add(&est); + return WriteContainer(cont, fNameOut); +} Index: trunk/MagicSoft/Mars/mjobs/MJOptimize.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimize.h (revision 6282) +++ trunk/MagicSoft/Mars/mjobs/MJOptimize.h (revision 6282) @@ -0,0 +1,191 @@ +#ifndef MARS_MJOptimize +#define MARS_MJOptimize + +#ifndef MARS_MJob +#include "MJob.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class TMinuit; + +class MAlphaFitter; + +class MEvtLoop; +class MParList; +class MFilter; +class MFilterList; +class MFitParameters; +class MFitParametersCalc; + +class MHMatrix; +class MGeomCam; +class MRead; +class MReadTree; + +class MJOptimize : public MJob +{ +public: + enum Optimizer_t + { + kMigrad, // Minimize by the method of Migrad + kSimplex, // Minimize by the method of Simplex + kMinimize, // Migrad + Simplex (if Migrad fails) + kMinos, // Minos error determination + kImprove, // Local minimum search + kSeek, // Minimize by the method of Monte Carlo + kNone // Skip optimization + }; + +private: + Int_t fDebug; // -1 no output, 0 MJOptimize output, 1 PrintStatistics output + Int_t fNumEvents; + + TList fRules; + TList fFilter; + TList fPreCuts; + + TList fNamesOn; + TList fNamesOff; + + TString fNameOut; + + void AddPoint(TList *l, Int_t idx, Float_t val) const; + TList *GetPlots() const; + + void AddRulesToMatrix(MHMatrix &m) const; + void SetupFilters(MFilterList &list, MFilter *filter=0) const; + Bool_t FillMatrix(MReadTree &read, MParList &l, Bool_t userules=kFALSE); + + MEvtLoop *fEvtLoop; //! + + Bool_t AddSequences(MRead &read, TList &list) const; + + // Minuit Interface + static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); + Double_t Fcn(const TArrayD &par, TMinuit *minuit=0); + + Int_t Minuit(TMinuit &minuit, const char *cmd) const; + Int_t Migrad(TMinuit &minuit) const { return Minuit(minuit, "MIGRAD"); } + Int_t Simplex(TMinuit &minuit) const { return Minuit(minuit, "SIMPLEX"); } + Int_t Minimize(TMinuit &minuit) const { return Minuit(minuit, "MINIMIZE"); } + Int_t Seek(TMinuit &minuit) const { return Minuit(minuit, "SEEK"); } + Int_t Improve(TMinuit &minuit) const { return Minuit(minuit, "IMPROVE"); } + Int_t Minos(TMinuit &minuit) const { return Minuit(minuit, "MINOS"); } + + TArrayD fParameters; //! + TArrayD fLimLo; //! + TArrayD fLimUp; //! + TArrayD fStep; //! + + Optimizer_t fType; + UInt_t fNumMaxCalls; + Float_t fTolerance; + + Bool_t Optimize(MEvtLoop &evtloop); + +public: + MJOptimize(); + + // I/O + void AddSequenceOn(const char *fname, const char *dir=""); + void AddSequenceOff(const char *fname, const char *dir=""); + + void ResetSequences(); + + // Interface for filter cuts + Int_t AddParameter(const char *rule); + void AddFilter(const char *rule); + void AddPreCut(const char *rule); + + // Steering of optimization + void SetNumEvents(UInt_t n); + void SetDebug(UInt_t n); + void SetNameOut(const char *name="") { fNameOut = name; } + void SetOptimizer(Optimizer_t o); + void SetNumMaxCalls(UInt_t num=0) { fNumMaxCalls=num; } + void SetTolerance(Float_t tol=0) { fTolerance=tol; } + + // Parameter access + void SetParameters(const TArrayD &par); + void SetParameter(Int_t idx, Double_t start, Double_t lo=0, Double_t up=0, Double_t step=-1) + { + if (fParameters.GetSize()<=idx) + { + fParameters.Set(idx+1); + fLimLo.Set(idx+1); + fLimUp.Set(idx+1); + fStep.Set(idx+1); + } + + fParameters[idx] = start; + fLimLo[idx] = lo; + fLimUp[idx] = up; + if (step<=0) + fStep[idx] = start==0 ? 0.1 : TMath::Abs(start*0.15); + else + fStep[idx] = step; + } + void FixParameter(Int_t idx, Double_t start) + { + if (fParameters.GetSize()<=idx) + { + fParameters.Set(idx+1); + fLimLo.Set(idx+1); + fLimUp.Set(idx+1); + fStep.Set(idx+1); + } + + fParameters[idx] = start; + fLimLo[idx] = 0; + fLimUp[idx] = 0; + fStep[idx] = 0; + } + + const TArrayD &GetParameters() const { return fParameters; } + + // Generalized optimizing routines + Bool_t Optimize(MParList &list); + + // Special optimizing routines + Bool_t Run(const char *fname, MFilter *filter, MAlphaFitter *fit=0); + Bool_t Run(const char *fname, MAlphaFitter *fit=0) + { + return Run(fname, 0, fit); + } + + Bool_t Run(MFilter *filter, MAlphaFitter *fit=0) + { + return Run(0, filter, fit); + } + Bool_t Run(MAlphaFitter *fit=0) + { + return Run(0, 0, fit); + } + + Bool_t RunOnOff(const char *fname, MFilter *filter, MAlphaFitter *fit=0, const char *tree="Events"); + Bool_t RunOnOff(const char *fname, MAlphaFitter *fit=0, const char *tree="Events") + { + return RunOnOff(fname, 0, fit, tree); + } + Bool_t RunEnergy(const char *fname, const char *rule); + + Bool_t RunOnOff(MFilter *filter, MAlphaFitter *fit=0, const char *tree="Events") + { + return RunOnOff(0, filter, fit, tree); + } + Bool_t RunOnOff(MAlphaFitter *fit=0, const char *tree="Events") + { + return RunOnOff(fit, tree); + } + Bool_t RunEnergy(const char *rule) + { + return RunEnergy(0, rule); + } + + ClassDef(MJOptimize, 0) // Class for optimization of the Supercuts +}; + +#endif Index: trunk/MagicSoft/Mars/mjobs/MSequences.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MSequences.cc (revision 6282) +++ trunk/MagicSoft/Mars/mjobs/MSequences.cc (revision 6282) @@ -0,0 +1,262 @@ +/* ======================================================================== *\ +! +! * +! * 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, 1/2005 +! +! Copyright: MAGIC Software Development, 2004-2005 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MSequences +// +// This class describes a collection of sequences. +// +// Such an input file looks like: +// +// crab.seq: +// --------- +// AnalysisNumber: 1 +// +// SequencesOn: 35222 +// SequencesOff: 36817 +// +// Sequence00035222.File: sequences/sequence035222.txt +// Sequence00036817.File: sequences/sequence036817.txt +// +// Sequence00035222.Dir: /data2/wuerzburg/Crab-Analyse/images/035222 +// Sequence00036817.Dir: /data2/wuerzburg/Crab-Analyse/images/036817 +// +// The analysis number is an artifical number used to name the output +// files automatically if the names are not overwritten in the corresponding +// programs. +// +// The sequence number are used to concatenate the filenames of the +// sequences using the file structure used in the datacenter. +// +// If you have different file names you can overwrite the default file names +// using Sequence%08d.File (make sure you have 8 digits!) +// +// In standard coditions (datacenter file system) paths are concatenated +// by using the information in the sequence files (date, etc). You can +// overwrite the directories in which the sequence-files (eg I-files) are +// stored using Sequence%08d.Dir (make sure you have 8 digits!) +// +// Resource file entries are case sensitive! +// +// MISSING (27/01/04): The default name and paths cannot be used yet, because +// they have to be defined soon. +// +///////////////////////////////////////////////////////////////////////////// +#include "MSequences.h" + +#include + +#include +#include +#include // TSystem::ExpandPath + +#include "MLog.h" +#include "MLogManip.h" + +#include "MRead.h" +#include "MDirIter.h" +#include "MSequence.h" + +ClassImp(MSequences); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Copy the run numbers from the TString runs into the TArrayI data +// +void MSequences::Split(TString &runs, TArrayI &data) const +{ + const TRegexp regexp("[0-9]+"); + + data.Set(0); + runs = runs.Strip(TString::kTrailing); + + while (!runs.IsNull()) + { + TString num = runs(regexp); + + const Int_t n = data.GetSize(); + data.Set(n+1); + data[n] = atoi(num.Data()); + + runs.Remove(0, runs.First(num)+num.Length()); + } +} + +void MSequences::ResolveSequences(TEnv &env, const TArrayI &num, TList &list) const +{ + for (int i=0; iAccessPathName(name, kFileExists)) + gLog << warn << "WARNING - Sequence file '" << name << "' doesn't exist." << endl; + + if (gSystem->AccessPathName(dir, kFileExists)) + gLog << warn << "WARNING - Directory '" << dir << "' doesn't exist." << endl; + + list.Add(new TNamed(name, dir)); + } +} + +// -------------------------------------------------------------------------- +// +// Read the file fname as setup file for the sequence. +// +MSequences::MSequences(const char *fname) +{ + fName = fname; + + const char *expname = gSystem->ExpandPathName(fname); + + fTitle = Form("Sequences contained in file %s", expname); + + TEnv env(expname); + delete [] expname; + + TString str; + + fNumAnalysis = env.GetValue("AnalysisNumber", -1); + + str = env.GetValue("SequencesOn", ""); + Split(str, fNumSequencesOn); + str = env.GetValue("SequencesOff", ""); + Split(str, fNumSequencesOff); + + + ResolveSequences(env, fNumSequencesOn, fSequencesOn); + ResolveSequences(env, fNumSequencesOff, fSequencesOff); + + + + //Print(); + /* + GetFileNames(env, fSequencesOn); + GetFileNames(env, fSequencesOff); + */ +} + +// -------------------------------------------------------------------------- +// +// Return '+' if both can be accessed, '-' otherwise. +// +void MSequences::PrintFile(const TObject &obj) +{ + const Bool_t access = !gSystem->AccessPathName(obj.GetName(), kFileExists) && !gSystem->AccessPathName(obj.GetTitle(), kFileExists) ? '+' : '-'; + gLog << " " << (access?"+":"-") << " " << obj.GetName() << " <" << obj.GetTitle() << ">" << endl; +} + +// -------------------------------------------------------------------------- +// +// Print the contents of the sequence +// +void MSequences::Print(Option_t *o) const +{ + gLog << all; + if (!IsValid()) + { + gLog << "Sequence: " << fName << " " << endl; + return; + } + gLog << "Analysis Number: " << fNumAnalysis << endl; + gLog << "Sequences On: "; + for (int i=0; i(&list)); + TObject *o=0; + while ((o=Next())) + { + MSequence seq(o->GetName()); + if (!seq.IsValid()) + { + gLog << warn << "WARNING - Sequence " << o->GetName() << " invalid!" << endl; + return kFALSE; + } + + const TString dir(o->GetTitle()); + seq.SetupDatRuns(files, dir.IsNull() ? 0 : dir.Data(), id, raw); + } + + return read.AddFiles(files)>0; +} + +Bool_t MSequences::AddFiles(MRead &read, char *id, Bool_t raw) const +{ + const Bool_t rc1 = AddFilesOff(read, id, raw); + const Bool_t rc2 = AddFilesOn(read, id, raw); + return rc1 && rc2; +} + +Bool_t MSequences::AddFilesOn(MRead &read, char *id, Bool_t raw) const +{ + return AddSequencesToList(fSequencesOn, read, id, raw); +} + +Bool_t MSequences::AddFilesOff(MRead &read, char *id, Bool_t raw) const +{ + return AddSequencesToList(fSequencesOff, read, id, raw); +} Index: trunk/MagicSoft/Mars/mjobs/MSequences.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MSequences.h (revision 6282) +++ trunk/MagicSoft/Mars/mjobs/MSequences.h (revision 6282) @@ -0,0 +1,55 @@ +#ifndef MARS_MSequences +#define MARS_MSequences + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayI +#include +#endif + +class MRead; + +class MSequences : public MParContainer +{ +private: + UInt_t fNumAnalysis; // Analysis number (artificial) + + TArrayI fNumSequencesOn; // number of on-sequences + TArrayI fNumSequencesOff; // number of off-sequences + + TList fSequencesOn; // list of names and paths of on-sequences + TList fSequencesOff; // list of names and paths of off-sequences + + void Split(TString &runs, TArrayI &data) const; + void ResolveSequences(TEnv &env, const TArrayI &num, TList &list) const; + static void PrintFile(const TObject &obj); + +public: + MSequences() : fNumAnalysis((UInt_t)-1) { } + MSequences(const char *fname); + + void Print(Option_t *o="") const; + + Bool_t IsValid() const { return fNumAnalysis!=(UInt_t)-1; } + + static Bool_t AddSequencesToList(const TList &list, MRead &read, char *id="I", Bool_t raw=kFALSE); + + Bool_t AddFiles(MRead &read, char *id="I", Bool_t raw=kFALSE) const; + Bool_t AddFilesOn(MRead &read, char *id="I", Bool_t raw=kFALSE) const; + Bool_t AddFilesOff(MRead &read, char *id="I", Bool_t raw=kFALSE) const; + + UInt_t GetNumSequencesOn() const { return fNumSequencesOn.GetSize(); } + UInt_t GetNumSequencesOff() const { return fNumSequencesOff.GetSize(); } + + Bool_t HasOffSequences() const { return GetNumSequencesOff()>0; } + + // Getter + UInt_t GetNumAnalysis() const { return fNumAnalysis; } + void SetNumAnalysis(UInt_t num) { fNumAnalysis=num; } + + ClassDef(MSequences, 0) +}; + +#endif