Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2662) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2663) @@ -48,4 +48,5 @@ #pragma link C++ class MCT1PadSchweizer+; #pragma link C++ class MCT1PadONOFF+; +#pragma link C++ class MPadONOFF+; #pragma link C++ class MCT1PointingCorrCalc+; @@ -60,4 +61,7 @@ #pragma link C++ class MCT1SupercutsCalc+; #pragma link C++ class MCT1FindSupercuts+; +#pragma link C++ class MSupercuts+; +#pragma link C++ class MSupercutsCalc+; +#pragma link C++ class MFindSupercuts+; #pragma link C++ class MMinuitInterface+; #pragma link C++ class MFiltercutsCalc+; @@ -68,4 +72,5 @@ #pragma link C++ class MCalibrationCam+; #pragma link C++ class MCalibrationCalc+; +#pragma link C++ class MPedestalWorkaround+; #pragma link C++ class MExtractedSignalCam+; Index: trunk/MagicSoft/Mars/manalysis/MFindSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MFindSupercuts.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MFindSupercuts.cc (revision 2663) @@ -0,0 +1,1173 @@ +/* ======================================================================== *\ +! +! * +! * 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, 6/2003 +! Author(s): Wolfgang Wittek, 7/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MFindSupercuts // +// // +// Class for otimizing the parameters of the supercuts // +// // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MFindSupercuts.h" + +#include // fabs + +#include +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MContinue.h" +#include "MSupercuts.h" +#include "MSupercutsCalc.h" +#include "MDataElement.h" +#include "MDataMember.h" + +#include "MEvtLoop.h" +#include "MFSelFinal.h" +#include "MF.h" +#include "MFEventSelector.h" +#include "MFEventSelector2.h" +#include "MFillH.h" +//#include "MGeomCamCT1Daniel.h" +#include "MFRandomSplit.h" +#include "MGeomCamMagic.h" +#include "MH3.h" +#include "MHSupercuts.h" +#include "MHFindSignificance.h" +#include "MHMatrix.h" +#include "MHOnSubtraction.h" + +#include "MLog.h" +#include "MLogManip.h" +#include "MMatrixLoop.h" +#include "MMinuitInterface.h" +#include "MParList.h" +#include "MProgressBar.h" +#include "MReadMarsFile.h" +#include "MReadTree.h" +#include "MTaskList.h" + + +ClassImp(MFindSupercuts); + +using namespace std; + + +//------------------------------------------------------------------------ +// +// fcnSupercuts +// +// - calculates the quantity to be minimized (using TMinuit) +// +// - the quantity to be minimized is (-1)*significance of the gamma signal +// in the alpha distribution (after cuts) +// +// - the parameters to be varied in the minimization are the cut parameters +// (par) +// +static void fcnSupercuts(Int_t &npar, Double_t *gin, Double_t &f, + Double_t *par, Int_t iflag) +{ + //cout << "entry fcnSupercuts" << endl; + + //------------------------------------------------------------- + // save pointer to the MINUIT object for optimizing the supercuts + // because it will be overwritten + // when fitting the alpha distribution in MHFindSignificance + TMinuit *savePointer = gMinuit; + //------------------------------------------------------------- + + + MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit(); + + MParList *plistfcn = (MParList*)evtloopfcn->GetParList(); + //MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList"); + + MSupercuts *super = (MSupercuts*)plistfcn->FindObject("MSupercuts"); + if (!super) + { + gLog << "fcnSupercuts : MSupercuts object '" << "MSupercuts" + << "' not found... aborting" << endl; + return; + } + + // + // transfer current parameter values to MSupercuts + // + // Attention : npar is the number of variable parameters + // not the total number of parameters + // + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); + + super->SetParameters(TArrayD(fNparx, par)); + + //$$$$$$$$$$$$$$$$$$$$$ + // for testing + //TArrayD checkparameters = super->GetParameters(); + //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", " + // << fNparx << endl; + //gLog << "fcnsupercuts : i, par, checkparameters =" << endl; + //for (Int_t i=0; iEventloop(); + + //tasklistfcn->PrintStatistics(0, kTRUE); + + MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3"); + if (!alpha) + return; + + TH1 &alphaHist = alpha->GetHist(); + alphaHist.SetName("alpha-fcnSupercuts"); + + //------------------------------------------- + // set Minuit pointer to zero in order not to destroy the TMinuit + // object for optimizing the supercuts + gMinuit = NULL; + + //================================================================= + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + const Double_t alphasig = 20.0; + const Double_t alphamin = 30.0; + const Double_t alphamax = 90.0; + const Int_t degree = 2; + + Bool_t drawpoly; + Bool_t fitgauss; + if (iflag == 3) + { + drawpoly = kTRUE; + fitgauss = kTRUE; + } + else + { + drawpoly = kFALSE; + fitgauss = kFALSE; + } + //drawpoly = kFALSE; + //fitgauss = kFALSE; + + const Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + const Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + + //================================================================= + + // reset gMinuit to the MINUIT object for optimizing the supercuts + gMinuit = savePointer; + //------------------------------------------- + + if (!rc) + { + gLog << "fcnSupercuts : FindSigma() failed" << endl; + f = 1.e10; + return; + } + + // plot some quantities during the optimization + MHSupercuts *plotsuper = (MHSupercuts*)plistfcn->FindObject("MHSupercuts"); + if (plotsuper) + plotsuper->Fill(&findsig); + + //------------------------ + // get significance + const Double_t significance = findsig.GetSignificance(); + f = significance>0 ? -significance : 0; + + + //------------------------ + // optimize signal/background ratio + //Double_t ratio = findsig.GetNbg()>0.0 ? + // findsig.GetNex()/findsig.GetNbg() : 0.0; + //f = -ratio; + + //------------------------------------------- + // final calculations + //if (iflag == 3) + //{ + // + //} + + //------------------------------------------------------------- +} + + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MFindSupercuts::MFindSupercuts(const char *name, const char *title) +: fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL) +{ + fName = name ? name : "MFindSupercuts"; + fTitle = title ? title : "Optimizer of the supercuts"; + + fUseOrigDistribution = kFALSE; + + //--------------------------- + // camera geometry is needed for conversion mm ==> degree + //fCam = new MGeomCamCT1Daniel; + fCam = new MGeomCamMagic; + + // matrices to contain the training/test samples + fMatrixTrain = new MHMatrix("MatrixTrain"); + fMatrixTest = new MHMatrix("MatrixTest"); + + // objects of MSupercutsCalc to which these matrices are attached + fCalcHadTrain = new MSupercutsCalc("SupercutsCalcTrain"); + fCalcHadTest = new MSupercutsCalc("SupercutsCalcTest"); + + // Define columns of matrices + fCalcHadTrain->InitMapping(fMatrixTrain); + fCalcHadTest->InitMapping(fMatrixTest); +} + +// -------------------------------------------------------------------------- +// +// Default destructor. +// +MFindSupercuts::~MFindSupercuts() +{ + delete fCam; + delete fMatrixTrain; + delete fMatrixTest; + delete fCalcHadTrain; + delete fCalcHadTest; +} + + +// -------------------------------------------------------------------------- +// +// +// +void MFindSupercuts::SetUseOrigDistribution(Bool_t b) +{ + fUseOrigDistribution = b; + + *fLog << "MFindSupercuts : when defining training and test matrices use the original distribution" + << endl; +} + +// -------------------------------------------------------------------------- +// +// Define the matrix 'fMatrixTrain' for the training sample +// +// alltogether 'howmanytrain' events are read from file 'nametrain'; +// the events are selected according to a target distribution 'hreftrain' +// +// +Bool_t MFindSupercuts::DefineTrainMatrix( + const TString &nametrain, MH3 &hreftrain, + const Int_t howmanytrain, const TString &filetrain) +{ + if (nametrain.IsNull() || howmanytrain <= 0) + return kFALSE; + + *fLog << "=============================================" << endl; + *fLog << "fill training matrix from file '" << nametrain + << "', select " << howmanytrain + << " events " << endl; + if (!fUseOrigDistribution) + { + *fLog << " according to a distribution given by the MH3 object '" + << hreftrain.GetName() << "'" << endl; + } + else + { + *fLog << " randomly" << endl; + } + + + MParList plist; + MTaskList tlist; + + MReadMarsFile read("Events", nametrain); + read.DisableAutoScheme(); + + MFEventSelector2 seltrain(hreftrain); + seltrain.SetNumMax(howmanytrain); + seltrain.SetName("selectTrain"); + if (fUseOrigDistribution) + seltrain.SetUseOrigDistribution(kTRUE); + + MFillH filltrain(fMatrixTrain); + filltrain.SetFilter(&seltrain); + filltrain.SetName("fillMatrixTrain"); + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTrain); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&seltrain); + tlist.AddToList(&filltrain); + + //****************************** + + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTrain"); + evtloop.SetProgressBar(&bar); + + if (!evtloop.Eventloop()) + return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + fMatrixTrain->Print("SizeCols"); + Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows(); + if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + { + *fLog << "MFindSupercuts::DefineTrainMatrix; no.of generated events (" + << howmanygenerated + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + *fLog << "training matrix was filled" << endl; + *fLog << "=============================================" << endl; + + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MFindSupercuts::DefineTrainMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Define the matrix for the test sample +// +// alltogether 'howmanytest' events are read from file 'nametest' +// +// the events are selected according to a target distribution 'hreftest' +// +// +Bool_t MFindSupercuts::DefineTestMatrix( + const TString &nametest, MH3 &hreftest, + const Int_t howmanytest, const TString &filetest) +{ + if (nametest.IsNull() || howmanytest<=0) + return kFALSE; + + *fLog << "=============================================" << endl; + *fLog << "fill test matrix from file '" << nametest + << "', select " << howmanytest + << " events " << endl; + if (!fUseOrigDistribution) + { + *fLog << " according to a distribution given by the MH3 object '" + << hreftest.GetName() << "'" << endl; + } + else + { + *fLog << " randomly" << endl; + } + + + MParList plist; + MTaskList tlist; + + MReadMarsFile read("Events", nametest); + read.DisableAutoScheme(); + + MFEventSelector2 seltest(hreftest); + seltest.SetNumMax(howmanytest); + seltest.SetName("selectTest"); + if (fUseOrigDistribution) + seltest.SetUseOrigDistribution(kTRUE); + + MFillH filltest(fMatrixTest); + filltest.SetFilter(&seltest); + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTest); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&seltest); + tlist.AddToList(&filltest); + + //****************************** + + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTest"); + evtloop.SetProgressBar(&bar); + + if (!evtloop.Eventloop()) + return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + fMatrixTest->Print("SizeCols"); + const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows(); + if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest)) + { + *fLog << "MFindSupercuts::DefineTestMatrix; no.of generated events (" + << howmanygenerated + << ") is incompatible with the no.of requested events (" + << howmanytest << ")" << endl; + } + + *fLog << "test matrix was filled" << endl; + *fLog << "=============================================" << endl; + + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MFindSupercuts::DefineTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Define the matrices for the training and test sample respectively +// +// +// +Bool_t MFindSupercuts::DefineTrainTestMatrix( + const TString &name, MH3 &href, + const Int_t howmanytrain, const Int_t howmanytest, + const TString &filetrain, const TString &filetest) +{ + *fLog << "=============================================" << endl; + *fLog << "fill training and test matrix from file '" << name + << "', select " << howmanytrain + << " training and " << howmanytest << " test events " << endl; + if (!fUseOrigDistribution) + { + *fLog << " according to a distribution given by the MH3 object '" + << href.GetName() << "'" << endl; + } + else + { + *fLog << " randomly" << endl; + } + + + MParList plist; + MTaskList tlist; + + MReadMarsFile read("Events", name); + read.DisableAutoScheme(); + + MFEventSelector2 selector(href); + selector.SetNumMax(howmanytrain+howmanytest); + selector.SetName("selectTrainTest"); + selector.SetInverted(); + if (fUseOrigDistribution) + selector.SetUseOrigDistribution(kTRUE); + + MContinue cont(&selector); + cont.SetName("ContTrainTest"); + + Double_t prob = ( (Double_t) howmanytrain ) + / ( (Double_t)(howmanytrain+howmanytest) ); + MFRandomSplit split(prob); + + MFillH filltrain(fMatrixTrain); + filltrain.SetFilter(&split); + filltrain.SetName("fillMatrixTrain"); + + + // consider this event as candidate for a test event + // only if event was not accepted as a training event + + MContinue conttrain(&split); + conttrain.SetName("ContTrain"); + + MFillH filltest(fMatrixTest); + filltest.SetName("fillMatrixTest"); + + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTrain); + plist.AddToList(fMatrixTest); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&cont); + + tlist.AddToList(&split); + tlist.AddToList(&filltrain); + tlist.AddToList(&conttrain); + + tlist.AddToList(&filltest); + + //****************************** + + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetName("EvtLoopMatrixTrainTest"); + evtloop.SetProgressBar(&bar); + + Int_t maxev = -1; + if (!evtloop.Eventloop(maxev)) + return kFALSE; + + tlist.PrintStatistics(0, kTRUE); + + fMatrixTrain->Print("SizeCols"); + const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows(); + if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + { + *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events (" + << generatedtrain + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + fMatrixTest->Print("SizeCols"); + const Int_t generatedtest = fMatrixTest->GetM().GetNrows(); + if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest)) + { + *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events (" + << generatedtest + << ") is incompatible with the no.of requested events (" + << howmanytest << ")" << endl; + } + + + *fLog << "training and test matrix were filled" << endl; + *fLog << "=============================================" << endl; + + + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MFindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MFindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Read training and test matrices from files +// +// + +Bool_t MFindSupercuts::ReadMatrix(const TString &filetrain, const TString &filetest) +{ + //-------------------------- + // read in training matrix + + TFile filetr(filetrain); + fMatrixTrain->Read("MatrixTrain"); + fMatrixTrain->Print("SizeCols"); + + *fLog << "MFindSupercuts::ReadMatrix; Training matrix was read in from file '" + << filetrain << "'" << endl; + filetr.Close(); + + + //-------------------------- + // read in test matrix + + TFile filete(filetest); + fMatrixTest->Read("MatrixTest"); + fMatrixTest->Print("SizeCols"); + + *fLog << "MFindSupercuts::ReadMatrix; Test matrix was read in from file '" + << filetest << "'" << endl; + filete.Close(); + + return kTRUE; +} + + +//------------------------------------------------------------------------ +// +// Steering program for optimizing the supercuts +// --------------------------------------------- +// +// the criterion for the 'optimum' is +// +// - a maximum significance of the gamma signal in the alpha plot, +// in which the supercuts have been applied +// +// The various steps are : +// +// - setup the event loop to be executed for each call to fcnSupercuts +// - call TMinuit to do the minimization of (-significance) : +// the fcnSupercuts function calculates the significance +// for the current cut values +// for this - the alpha plot is produced in the event loop, +// in which the cuts have been applied +// - the significance of the gamma signal in the alpha plot +// is calculated +// +// Needed as input : (to be set by the Set functions) +// +// - fFilenameParam name of file to which optimum values of the +// parameters are written +// +// - fHadronnessName name of container where MSupercutsCalc will +// put the supercuts hadronness +// +// - for the minimization, the starting values of the parameters are taken +// - from file parSCinit (if it is != "") +// - or from the arrays params and/or steps +// +//---------------------------------------------------------------------- +Bool_t MFindSupercuts::FindParams(TString parSCinit, + TArrayD ¶ms, TArrayD &steps) +{ + // Setup the event loop which will be executed in the + // fcnSupercuts function of MINUIT + // + // parSCinit is the name of the file containing the initial values + // of the parameters; + // if parSCinit = "" 'params' and 'steps' are taken as initial values + + *fLog << "=============================================" << endl; + *fLog << "Setup event loop for fcnSupercuts" << endl; + + if (fHadronnessName.IsNull()) + { + *fLog << "MFindSupercuts::FindParams; hadronness name is not defined... aborting" + << endl; + return kFALSE; + } + + if (fMatrixTrain == NULL) + { + *fLog << "MFindSupercuts::FindParams; training matrix is not defined... aborting" + << endl; + return kFALSE; + } + + if (fMatrixTrain->GetM().GetNrows() <= 0) + { + *fLog << "MFindSupercuts::FindParams; training matrix has no entries" + << endl; + return kFALSE; + } + + MParList parlistfcn; + MTaskList tasklistfcn; + + // loop over rows of matrix + MMatrixLoop loop(fMatrixTrain); + + //-------------------------------- + // create container for the supercut parameters + // and set them to their initial values + MSupercuts super; + + // take initial values from file parSCinit + if (parSCinit != "") + { + TFile inparam(parSCinit); + super.Read("MSupercuts"); + inparam.Close(); + *fLog << "MFindSupercuts::FindParams; initial values of parameters are taken from file " + << parSCinit << endl; + } + + // take initial values from 'params' and/or 'steps' + else if (params.GetSize() != 0 || steps.GetSize() != 0 ) + { + if (params.GetSize() != 0) + { + *fLog << "MFindSupercuts::FindParams; initial values of parameters are taken from 'params'" + << endl; + super.SetParameters(params); + } + if (steps.GetSize() != 0) + { + *fLog << "MFindSupercuts::FindParams; initial step sizes are taken from 'steps'" + << endl; + super.SetStepsizes(steps); + } + } + else + { + *fLog << "MFindSupercuts::FindParams; initial values and step sizes are taken from the MSupercuts constructor" + << endl; + } + + + //-------------------------------- + // calculate supercuts hadronness + fCalcHadTrain->SetHadronnessName(fHadronnessName); + + // apply the supercuts + MF scfilter(fHadronnessName+".fHadronness>0.5"); + MContinue supercuts(&scfilter); + + // plot |alpha| + const TString mh3Name = "AlphaFcn"; + MBinning binsalpha("Binning"+mh3Name); + binsalpha.SetEdges(54, -12.0, 96.0); + + *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl; + + // |alpha| is assumed to be in column 7 of the matrix + MH3 alpha("MatrixTrain[7]"); + alpha.SetName(mh3Name); + + MFillH fillalpha(&alpha); + + // book histograms to be filled during the optimization + // ! not in the event loop ! + MHSupercuts plotsuper; + plotsuper.SetupFill(&parlistfcn); + + //****************************** + // entries in MParList (extension of old MParList) + + parlistfcn.AddToList(&tasklistfcn); + parlistfcn.AddToList(&super); + parlistfcn.AddToList(fCam); + parlistfcn.AddToList(fMatrixTrain); + + parlistfcn.AddToList(&binsalpha); + parlistfcn.AddToList(&alpha); + + parlistfcn.AddToList(&plotsuper); + + //****************************** + // entries in MTaskList + + tasklistfcn.AddToList(&loop); + tasklistfcn.AddToList(fCalcHadTrain); + tasklistfcn.AddToList(&supercuts); + tasklistfcn.AddToList(&fillalpha); + + + //****************************** + + MEvtLoop evtloopfcn("EvtLoopFCN"); + evtloopfcn.SetParList(&parlistfcn); + *fLog << "Event loop for fcnSupercuts has been setup" << endl; + + // address of evtloopfcn is used in CallMinuit() + + + //----------------------------------------------------------------------- + // + //---------- Start of minimization part -------------------- + // + // Do the minimization with MINUIT + // + // Be careful: This is not thread safe + // + *fLog << "========================================================" << endl; + *fLog << "Start minimization for supercuts" << endl; + + + // ------------------------------------------- + // prepare call to MINUIT + // + + // get initial values of parameters + fVinit = super.GetParameters(); + fStep = super.GetStepsizes(); + + TString name[fVinit.GetSize()]; + fStep.Set(fVinit.GetSize()); + fLimlo.Set(fVinit.GetSize()); + fLimup.Set(fVinit.GetSize()); + fFix.Set(fVinit.GetSize()); + + fNpar = fVinit.GetSize(); + + for (UInt_t i=0; i= 48) + // { + // fStep[i] = 0.0; + // fFix[i] = 1; + // } + //} + + + // ------------------------------------------- + // call MINUIT + + TStopwatch clock; + clock.Start(); + + *fLog << "before calling CallMinuit" << endl; + + MMinuitInterface inter; + Bool_t rc = inter.CallMinuit(fcnSupercuts, name, + fVinit, fStep, fLimlo, fLimup, fFix, + &evtloopfcn, "SIMPLEX", kFALSE); + + *fLog << "after calling CallMinuit" << endl; + + *fLog << "Time spent for the minimization in MINUIT : " << endl;; + clock.Stop(); + clock.Print(); + + plotsuper.DrawClone(); + + if (!rc) + return kFALSE; + + *fLog << "Minimization for supercuts finished" << endl; + *fLog << "========================================================" << endl; + + + // ----------------------------------------------------------------- + // in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put + // into MSupercuts + + // write optimum parameter values onto file filenameParam + + TFile outparam(fFilenameParam, "RECREATE"); + super.Write(); + outparam.Close(); + + *fLog << "Optimum parameter values for supercuts were written onto file '" + << fFilenameParam << "' :" << endl; + + const TArrayD &check = super.GetParameters(); + for (Int_t i=0; iGetM().GetNrows() <= 0) + { + *fLog << "MFindSupercuts::TestParams; test matrix has no entries" + << endl; + return kFALSE; + } + + // ------------- TEST the supercuts ------------------------------ + // + *fLog << "Test the supercuts on the test sample" << endl; + + // ----------------------------------------------------------------- + // read optimum parameter values from file filenameParam + // into array 'supercutsPar' + + TFile inparam(fFilenameParam); + MSupercuts scin; + scin.Read("MSupercuts"); + inparam.Close(); + + *fLog << "Optimum parameter values for supercuts were read from file '"; + *fLog << fFilenameParam << "' :" << endl; + + const TArrayD &supercutsPar = scin.GetParameters(); + for (Int_t i=0; iSetHadronnessName(fHadronnessName); + + + //------------------------------------------- + + MMatrixLoop loop(fMatrixTest); + + // plot alpha before applying the supercuts + const TString mh3NameB = "AlphaBefore"; + MBinning binsalphabef("Binning"+mh3NameB); + binsalphabef.SetEdges(54, -12.0, 96.0); + + // |alpha| is assumed to be in column 7 of the matrix + MH3 alphabefore("MatrixTest[7]"); + alphabefore.SetName(mh3NameB); + + TH1 &alphahistb = alphabefore.GetHist(); + alphahistb.SetName("alphaBefore-TestParams"); + + MFillH fillalphabefore(&alphabefore); + fillalphabefore.SetName("FillAlphaBefore"); + + // apply the supercuts + MF scfilter(fHadronnessName+".fHadronness>0.5"); + MContinue applysupercuts(&scfilter); + + // plot alpha after applying the supercuts + const TString mh3NameA = "AlphaAfter"; + MBinning binsalphaaft("Binning"+mh3NameA); + binsalphaaft.SetEdges(54, -12.0, 96.0); + + MH3 alphaafter("MatrixTest[7]"); + alphaafter.SetName(mh3NameA); + + TH1 &alphahista = alphaafter.GetHist(); + alphahista.SetName("alphaAfter-TestParams"); + + MFillH fillalphaafter(&alphaafter); + fillalphaafter.SetName("FillAlphaAfter"); + + //****************************** + // entries in MParList + + parlist2.AddToList(&tasklist2); + + parlist2.AddToList(&supercuts); + + parlist2.AddToList(fCam); + parlist2.AddToList(fMatrixTest); + + parlist2.AddToList(&binsalphabef); + parlist2.AddToList(&alphabefore); + + parlist2.AddToList(&binsalphaaft); + parlist2.AddToList(&alphaafter); + + //****************************** + // entries in MTaskList + + tasklist2.AddToList(&loop); + tasklist2.AddToList(&fillalphabefore); + + tasklist2.AddToList(fCalcHadTest); + tasklist2.AddToList(&applysupercuts); + + tasklist2.AddToList(&fillalphaafter); + + //****************************** + + MProgressBar bar2; + MEvtLoop evtloop2; + evtloop2.SetParList(&parlist2); + evtloop2.SetName("EvtLoopTestParams"); + evtloop2.SetProgressBar(&bar2); + Int_t maxevents2 = -1; + if (!evtloop2.Eventloop(maxevents2)) + return kFALSE; + + tasklist2.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // draw the alpha plots + + MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3"); + TH1 &alphaHist1 = alphaBefore->GetHist(); + alphaHist1.SetXTitle("|\\alpha| [\\circ]"); + + MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3"); + TH1 &alphaHist2 = alphaAfter->GetHist(); + alphaHist2.SetXTitle("|\\alpha| [\\circ]"); + alphaHist2.SetName("alpha-TestParams"); + + TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300); + c->Divide(2,1); + + gROOT->SetSelectedPad(NULL); + + c->cd(1); + alphaHist1.DrawCopy(); + + c->cd(2); + alphaHist2.DrawCopy(); + + + + //------------------------------------------- + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + const Double_t alphasig = 20.0; + const Double_t alphamin = 30.0; + const Double_t alphamax = 90.0; + const Int_t degree = 2; + const Bool_t drawpoly = kTRUE; + const Bool_t fitgauss = kTRUE; + const Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + + const Double_t significance = findsig.GetSignificance(); + const Double_t alphasi = findsig.GetAlphasi(); + + *fLog << "Significance of gamma signal after supercuts in test sample : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + //------------------------------------------- + + + *fLog << "Test of supercuts part finished" << endl; + *fLog << "======================================================" << endl; + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/manalysis/MFindSupercuts.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MFindSupercuts.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MFindSupercuts.h (revision 2663) @@ -0,0 +1,134 @@ +#ifndef MARS_MFindSupercuts +#define MARS_MFindSupercuts + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif +#ifndef ROOT_TArrayI +#include +#endif + +class MFilter; +class MEvtLoop; +class MH3; +class MSupercutsCalc; +class MGeomCam; +class MHMatrix; +/* +#include "MFilter.h" +#include "MEvtLoop.h" +#include "MH3.h" +#include "MSupercutsCalc.h" +#include "MGeomCam.h" +#include "MHMatrix.h" +*/ + +class MFindSupercuts : public MParContainer +{ +private: + + TString fFilenameTrain; + TString fFilenameTest; + + Int_t fHowManyTrain; + Int_t fHowManyTest; + + Bool_t fUseOrigDistribution; + + TString fFilenameParam; + + TString fHadronnessName; + + MSupercutsCalc *fCalcHadTrain; + MSupercutsCalc *fCalcHadTest; + + MHMatrix *fMatrixTrain; + MHMatrix *fMatrixTest; + MGeomCam *fCam; + + MEvtLoop *fObjectFit; + + MFilter *fMatrixFilter; + + // to comunicate with MINUIT ----------------- + // attention : dimensions must agree with those in + // MMinuitInterface::CallMinuit() + //char fParName [80][100]; + TArrayD fVinit; + TArrayD fStep; + TArrayD fLimlo; + TArrayD fLimup; + TArrayI fFix; + + UInt_t fNpar; + + TString fMethod; + + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + Int_t fErrMinimize; + //-------------------------------------------- + + +public: + MFindSupercuts(const char *name=NULL, const char *title=NULL); + ~MFindSupercuts(); + + void SetFilenameTraining(const TString &name, const Int_t howmany) + {fFilenameTrain = name; fHowManyTrain = howmany; } + + void SetFilenameTest(const TString &name, const Int_t howmany) + {fFilenameTest = name; fHowManyTest = howmany; } + + void SetFilenameParam(const TString &name) {fFilenameParam = name;} + void SetHadronnessName(const TString &name) {fHadronnessName = name;} + + void SetMatrixFilter(MFilter *filter) {fMatrixFilter = filter;} + + Bool_t DefineTrainMatrix(const TString &name, MH3 &href, + const Int_t howmany, const TString &filetrain); + + Bool_t DefineTestMatrix(const TString &name, MH3 &href, + const Int_t howmany, const TString &filetest); + + Bool_t DefineTrainTestMatrix(const TString &name, MH3 &href, + const Int_t howmanytrain, const Int_t howmanytest, + const TString &filetrain, const TString &filetest); + + void SetUseOrigDistribution(Bool_t b); + + MHMatrix *GetMatrixTrain() { return fMatrixTrain; } + MHMatrix *GetMatrixTest() { return fMatrixTest; } + + Bool_t ReadMatrix( const TString &filetrain, const TString &filetest); + + Bool_t FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps); + Bool_t TestParams(); + + ClassDef(MFindSupercuts, 1) // Class for optimization of the Supercuts +}; + +#endif + + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2663) @@ -0,0 +1,1874 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek, 06/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MPadONOFF +// +// This task applies padding such that for a given pixel and for a given +// Theta bin the resulting distribution of the pedestal sigma is identical +// to the distributions given by fHSigmaPixTheta and fHDiffPixTheta. +// +// The number of photons, its error and the pedestal sigmas are altered. +// On average, the number of photons added is zero. +// +// The formulas used can be found in Thomas Schweizer's Thesis, +// Section 2.2.1 +// +// There are 2 options for the padding : +// +// 1) fPadFlag = 1 : +// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta +// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using +// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 +// (fHDiffPixTheta). +// +// This is the preferred option as it takes into account the +// correlations between the Sigma of a pixel and Sigmabar. +// +// THE FOLLOWING OPTION HAS NOT BEEN TESTED : +// 2) fPadFlag = 2 : +// Generate a pedestal sigma for each pixel using the 3D-histogram +// Theta, pixel no., Sigma (fHSigmaPixTheta). +// +// +// The padding has to be done before the image cleaning because the +// image cleaning depends on the pedestal sigmas. +// +// For random numbers gRandom is used. +// +///////////////////////////////////////////////////////////////////////////// +#include "MPadONOFF.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MSigmabar.h" +#include "MMcEvt.hxx" +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MGeomCam.h" + +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MBlindPixels.h" + +#include "MRead.h" +#include "MFilterList.h" +#include "MTaskList.h" +#include "MBlindPixelCalc.h" +#include "MHBlindPixels.h" +#include "MFillH.h" +#include "MHSigmaTheta.h" +#include "MEvtLoop.h" + +ClassImp(MPadONOFF); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MPadONOFF::MPadONOFF(const char *name, const char *title) +{ + fName = name ? name : "MPadONOFF"; + fTitle = title ? title : "Task for the ON-OFF padding"; + + fPadFlag = 1; + *fLog << all << "fPadFlag = " << fPadFlag << endl; + + fType = ""; + + fHSigmaTheta = NULL; + fHSigmaThetaON = NULL; + fHSigmaThetaOFF = NULL; + + fHSigmaPixTheta = NULL; + fHSigmaPixThetaON = NULL; + fHSigmaPixThetaOFF = NULL; + + fHDiffPixTheta = NULL; + fHDiffPixThetaON = NULL; + fHDiffPixThetaOFF = NULL; + + fHgON = NULL; + fHgOFF = NULL; + + fHBlindPixIdTheta = NULL; + fHBlindPixNTheta = NULL; + + fHSigmaPedestal = NULL; + fHPhotons = NULL; + fHNSB = NULL; +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MPadONOFF::~MPadONOFF() +{ + if (fHBlindPixIdTheta != NULL) delete fHBlindPixIdTheta; + if (fHBlindPixNTheta != NULL) delete fHBlindPixNTheta; + + if (fHSigmaTheta != NULL) delete fHSigmaTheta; + if (fHSigmaPixTheta != NULL) delete fHSigmaPixTheta; + if (fHDiffPixTheta != NULL) delete fHDiffPixTheta; + + if (fHSigmaPedestal != NULL) delete fHSigmaPedestal; + if (fHPhotons != NULL) delete fHPhotons; + if (fHNSB != NULL) delete fHNSB; + + if (fInfile != NULL) delete fInfile; +} + +// -------------------------------------------------------------------------- +// +// Merge the distributions of ON and OFF data +// +// fHSigmaTheta 2D-histogram (Theta, sigmabar) +// fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma) +// fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) +// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id) +// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels ) +// +// and define the target distributions for the padding +// +Bool_t MPadONOFF::MergeHistograms(TH2D *sigthon, TH2D *sigthoff, + TH3D *sigpixthon, TH3D *sigpixthoff, + TH3D *diffpixthon, TH3D *diffpixthoff, + TH2D *blindidthon, TH2D *blindidthoff, + TH2D *blindnthon, TH2D *blindnthoff) +{ + *fLog << "----------------------------------------------------------------------------------" << endl; + *fLog << all << "Merge the ON and OFF histograms to obtain the target distributions for the padding" + << endl; + + //------------------------------------------------------------- + // merge the distributions of ON and OFF events + // to obtain the target distribution fHSigmaTheta + // + + Double_t eps = 1.e-10; + + fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); + fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); + + // get binning for fHgON and fHgOFF from sigthon + // binning in Theta + TAxis *ax = sigthon->GetXaxis(); + Int_t nbinstheta = ax->GetNbins(); + TArrayD edgesx; + edgesx.Set(nbinstheta+1); + for (Int_t i=0; i<=nbinstheta; i++) + { + edgesx[i] = ax->GetBinLowEdge(i+1); + *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] + << endl; + } + MBinning binth; + binth.SetEdges(edgesx); + + // binning in sigmabar + TAxis *ay = sigthon->GetYaxis(); + Int_t nbinssig = ay->GetNbins(); + TArrayD edgesy; + edgesy.Set(nbinssig+1); + for (Int_t i=0; i<=nbinssig; i++) + { + edgesy[i] = ay->GetBinLowEdge(i+1); + *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] + << endl; + } + MBinning binsg; + binsg.SetEdges(edgesy); + + + fHgON = new TH3D; + MH::SetBinning(fHgON, &binth, &binsg, &binsg); + fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); + + fHgOFF = new TH3D; + MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); + fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); + + //............ loop over Theta bins ........................... + + // hista is the normalized 1D histogram of sigmabar for ON data + // histb is the normalized 1D histogram of sigmabar for OFF data + // histc is the difference ON-OFF + + // at the beginning, histap is a copy of hista + // at the end, it will be the 1D histogram for ON data after padding + + // at the beginning, histbp is a copy of histb + // at the end, it will be the 1D histogram for OFF data after padding + + // at the beginning, histcp is a copy of histc + // at the end, it should be zero + + *fLog << all << "MergeHistograms; bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; + for (Int_t l=1; l<=nbinstheta; l++) + { + TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); + Stat_t suma = hista->Integral(); + hista->Scale(1./suma); + + TH1D *histap = new TH1D( (const TH1D&)*hista ); + + TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, ""); + Stat_t sumb = histb->Integral(); + histb->Scale(1./sumb); + + TH1D *histbp = new TH1D( (const TH1D&)*histb ); + + TH1D *histc = new TH1D( (const TH1D&)*hista ); + histc->Add(histb, -1.0); + + TH1D *histcp = new TH1D( (const TH1D&)*histc ); + + + // calculate matrix g + + // fHgON[k][j] tells which fraction of ON events should be padded + // from sigma_k to sigma_j + + + // fHgOFF[k][j] tells which fraction of OFF events should be padded + // from sigma_k to sigma_j + + //-------- start j loop ------------------------------------------------ + // loop over bins in histc, starting from the end + Double_t v, s, w, t, x, u, a, b, c, arg; + + for (Int_t j=nbinssig; j >= 1; j--) + { + v = histcp->GetBinContent(j); + if ( fabs(v) < eps ) continue; + if (v >= 0.0) + s = 1.0; + else + s = -1.0; + + //................ start k loop ...................................... + // look for a bin k which may compensate the content of bin j + for (Int_t k=j-1; k >= 1; k--) + { + w = histcp->GetBinContent(k); + if ( fabs(w) < eps ) continue; + if (w >= 0.0) + t = 1.0; + else + t = -1.0; + + if (s==t) continue; + + x = v + w; + if (x >= 0.0) + u = 1.0; + else + u = -1.0; + + if (u == s) + { + arg = -w; + + if (arg <=0.0) + { + fHgON->SetBinContent (l, k, j, -arg); + fHgOFF->SetBinContent(l, k, j, 0.0); + } + else + { + fHgON->SetBinContent (l, k, j, 0.0); + fHgOFF->SetBinContent(l, k, j, arg); + } + + + *fLog << "l, k, j, arg = " << l << ", " << k << ", " << j + << ", " << arg << endl; + + // g(k-1, j-1) = arg; + //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", " + // << arg << endl; + + //...................................... + // this is for checking the procedure + if (arg < 0.0) + { + a = histap->GetBinContent(k); + histap->SetBinContent(k, a+arg); + a = histap->GetBinContent(j); + histap->SetBinContent(j, a-arg); + } + else + { + b = histbp->GetBinContent(k); + histbp->SetBinContent(k, b-arg); + b = histbp->GetBinContent(j); + histbp->SetBinContent(j, b+arg); + } + //...................................... + + histcp->SetBinContent(k, 0.0); + histcp->SetBinContent(j, x); + + //...................................... + // redefine v + v = histcp->GetBinContent(j); + if ( fabs(v) < eps ) break; + if (v >= 0.0) + s = 1.0; + else + s = -1.0; + //...................................... + + continue; + } + + arg = v; + + if (arg <=0.0) + { + fHgON->SetBinContent (l, k, j, -arg); + fHgOFF->SetBinContent(l, k, j, 0.0); + } + else + { + fHgON->SetBinContent (l, k, j, 0.0); + fHgOFF->SetBinContent(l, k, j, arg); + } + + *fLog << "l, k, j, arg = " << l << ", " << k << ", " << j + << ", " << arg << endl; + + //g(k-1, j-1) = arg; + //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", " + // << arg << endl; + + //...................................... + // this is for checking the procedure + if (arg < 0.0) + { + a = histap->GetBinContent(k); + histap->SetBinContent(k, a+arg); + a = histap->GetBinContent(j); + histap->SetBinContent(j, a-arg); + } + else + { + b = histbp->GetBinContent(k); + + histbp->SetBinContent(k, b-arg); + b = histbp->GetBinContent(j); + histbp->SetBinContent(j, b+arg); + } + //...................................... + + histcp->SetBinContent(k, x); + histcp->SetBinContent(j, 0.0); + + break; + } + //................ end k loop ...................................... + } + //-------- end j loop ------------------------------------------------ + + // check results for this Theta bin + for (Int_t j=1; j<=nbinssig; j++) + { + a = histap->GetBinContent(j); + b = histbp->GetBinContent(j); + c = histcp->GetBinContent(j); + + if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps ) + *fLog << err << "Read; inconsistency in results; a, b, c = " + << a << ", " << b << ", " << c << endl; + } + + + // fill target distribution SigmaTheta + // for this Theta bin + // + for (Int_t j=1; j<=nbinssig; j++) + { + a = histap->GetBinContent(j); + fHSigmaTheta->SetBinContent(l, j, a); + } + + delete hista; + delete histb; + delete histc; + + delete histap; + delete histbp; + delete histcp; + } + //............ end of loop over Theta bins .................... + + + // target distributions for MC + // SigmaPixTheta and DiffPixTheta + // BlindPixIdTheta and BlindPixNTheta + // are calculated as averages of the ON and OFF distributions + + // Note : in all histograms the pedestal RMS is assumed to be given in units of "number of photons" + + fHSigmaThetaON = sigthon; + fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (target)"); + + fHSigmaThetaOFF = sigthoff; + fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (target)"); + + + fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); + fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); + + fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); + fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); + + fHSigmaPixTheta = new TH3D( (TH3D&)*sigpixthon ); + fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); + + fHSigmaPixThetaON = sigpixthon; + fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (target)"); + + fHSigmaPixThetaOFF = sigpixthoff; + fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (target)"); + + fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); + fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); + + fHDiffPixThetaON = diffpixthon; + fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (target)"); + + fHDiffPixThetaOFF = diffpixthoff; + fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (target)"); + + + for (Int_t j=1; j<=nbinstheta; j++) + { + // fraction of ON/OFF events to be padded to a sigmabar + // of the OFF/ON events : fracON, fracOFF + + Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + + TAxis *ax; + TAxis *ay; + TAxis *az; + + Double_t entON; + Double_t entOFF; + + Double_t normON; + Double_t normOFF; + + // set ranges for 2D-projections of 3D-histograms + ax = sigpixthon->GetXaxis(); + ax->SetRange(j, j); + ax = sigpixthoff->GetXaxis(); + ax->SetRange(j, j); + + ax = diffpixthon->GetXaxis(); + ax->SetRange(j, j); + ax = diffpixthoff->GetXaxis(); + ax->SetRange(j, j); + + TH2D *hist; + TH2D *histOFF; + + TH1D *hist1; + TH1D *hist1OFF; + + // weights for ON and OFF distrubtions when + // calculating the weighted average + Double_t facON = 0.5 * (1. - fracON + fracOFF); + Double_t facOFF = 0.5 * (1. + fracON - fracOFF); + + *fLog << all << fracON << ", " << fracOFF << ", " + << facON << ", " << facOFF << endl; + //------------------------------------------- + ay = blindnthon->GetYaxis(); + Int_t nbinsn = ay->GetNbins(); + + hist1 = (TH1D*)blindnthon->ProjectionY ("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = (TH1D*)blindnthoff->ProjectionY("", j, j, ""); + + // number of events in this theta bin + entON = hist1->Integral(); + entOFF = hist1OFF->Integral(); + + *fLog << all << "BlindPixNTheta : j, entON, entOFF = " << j << ", " + << entON << ", " << entOFF << endl; + + normON = entON<=0.0 ? 0.0 : 1.0/entON; + normOFF = entOFF<=0.0 ? 0.0 : 1.0/entOFF; + + hist1->Scale(normON); + hist1OFF->Scale(normOFF); + + // weighted average of ON and OFF distributions + hist1->Scale(facON); + hist1->Add(hist1OFF, facOFF); + + for (Int_t k=1; k<=nbinsn; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixNTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + + //------------------------------------------- + ay = blindidthon->GetYaxis(); + Int_t nbinsid = ay->GetNbins(); + + hist1 = (TH1D*)blindidthon->ProjectionY ("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = (TH1D*)blindidthoff->ProjectionY("", j, j, ""); + + hist1->Scale(normON); + hist1OFF->Scale(normOFF); + + // weighted average of ON and OFF distributions + hist1->Scale(facON); + hist1->Add(hist1OFF, facOFF); + + for (Int_t k=1; k<=nbinsid; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixIdTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + + //------------------------------------------- + ay = sigpixthon->GetYaxis(); + Int_t nbinspix = ay->GetNbins(); + + az = sigpixthon->GetZaxis(); + Int_t nbinssigma = az->GetNbins(); + + hist = (TH2D*)sigpixthon->Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)sigpixthoff->Project3D("zy"); + + hist->Scale(normON); + histOFF->Scale(normOFF); + + // weighted average of ON and OFF distributions + + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspix; k++) + for (Int_t l=1; l<=nbinssigma; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHSigmaPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + //------------------------------------------- + ay = diffpixthon->GetYaxis(); + Int_t nbinspixel = ay->GetNbins(); + + az = diffpixthon->GetZaxis(); + Int_t nbinsdiff = az->GetNbins(); + + hist = (TH2D*)diffpixthon->Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)diffpixthoff->Project3D("zy"); + + hist->Scale(normON); + histOFF->Scale(normOFF); + + // weighted average of ON and OFF distributions + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinsdiff; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHDiffPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + //------------------------------------------- + } + + + *fLog << all << "The target distributions for the padding have been created" + << endl; + *fLog << all << "----------------------------------------------------------" + << endl; + //-------------------------------------------- + + fHSigmaTheta->SetDirectory(NULL); + fHSigmaThetaON->SetDirectory(NULL); + fHSigmaThetaOFF->SetDirectory(NULL); + + fHSigmaPixTheta->SetDirectory(NULL); + fHSigmaPixThetaON->SetDirectory(NULL); + fHSigmaPixThetaOFF->SetDirectory(NULL); + + fHDiffPixTheta->SetDirectory(NULL); + fHDiffPixThetaON->SetDirectory(NULL); + fHDiffPixThetaOFF->SetDirectory(NULL); + + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetDirectory(NULL); + + + fHgON->SetDirectory(NULL); + fHgOFF->SetDirectory(NULL); + + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Read target distributions from a file +// +// +Bool_t MPadONOFF::ReadTargetDist(const char* namefilein) +{ + *fLog << all << "Read padding histograms from file " << namefilein << endl; + + fInfile = new TFile(namefilein); + fInfile->ls(); + + //------------------------------------ + + fHSigmaTheta = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabar"); + if (!fHSigmaTheta) + { + *fLog << all << "Object '2D-ThetaSigmabar' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaSigmabar' was read in" << endl; + + fHSigmaThetaON = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON"); + if (!fHSigmaThetaON) + { + *fLog << all << "Object '2D-ThetaSigmabarON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaSigmabarON' was read in" << endl; + + fHSigmaThetaOFF = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF"); + if (!fHSigmaThetaOFF) + { + *fLog << "Object '2D-ThetaSigmabarOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaSigmabarOFF' was read in" << endl; + + fHSigmaPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); + if (!fHSigmaPixTheta) + { + *fLog << all << "Object '3D-ThetaPixSigma' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-ThetaPixSigma' was read in" << endl; + + fHSigmaPixThetaON = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON"); + if (!fHSigmaPixThetaON) + { + *fLog << all << "Object '3D-ThetaPixSigmaON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-ThetaPixSigmaON' was read in" << endl; + + fHSigmaPixThetaOFF = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF"); + if (!fHSigmaPixThetaOFF) + { + *fLog << all << "Object '3D-ThetaPixSigmaOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << "Object '3D-ThetaPixSigmaOFF' was read in" << endl; + + fHDiffPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); + if (!fHDiffPixTheta) + { + *fLog << all << "Object '3D-ThetaPixDiff' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-ThetaPixDiff' was read in" << endl; + + fHDiffPixThetaON = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON"); + if (!fHDiffPixThetaON) + { + *fLog << all << "Object '3D-ThetaPixDiffON' not found on root file" + << endl; + return kFALSE; + } + *fLog << "Object '3D-ThetaPixDiffON' was read in" << endl; + + fHDiffPixThetaOFF = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF"); + if (!fHDiffPixThetaOFF) + { + *fLog << all << "Object '3D-ThetaPixDiffOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-ThetaPixDiffOFF' was read in" << endl; + + fHgON = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixON"); + if (!fHgON) + { + *fLog << all << "Object '3D-PaddingMatrixON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-PaddingMatrixON' was read in" << endl; + + fHgOFF = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF"); + if (!fHgOFF) + { + *fLog << all << "Object '3D-PaddingMatrixOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '3D-PaddingMatrixOFF' was read in" << endl; + + + fHBlindPixIdTheta = + (TH2D*) gROOT->FindObject("2D-ThetaBlindId"); + if (!fHBlindPixIdTheta) + { + *fLog << all << "Object '2D-ThetaBlindId' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaBlindId' was read in" << endl; + + fHBlindPixNTheta = + (TH2D*) gROOT->FindObject("2D-ThetaBlindN"); + if (!fHBlindPixNTheta) + { + *fLog << all << "Object '2D-ThetaBlindN' not found on root file" + << endl; + return kFALSE; + } + *fLog << all << "Object '2D-ThetaBlindN' was read in" << endl; + + + //------------------------------------ + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Write target distributions onto a file +// +// +Bool_t MPadONOFF::WriteTargetDist(const char* namefileout) +{ + *fLog << all << "Write padding histograms onto file " << namefileout << endl; + + TFile outfile(namefileout, "RECREATE"); + + fHBlindPixNTheta->Write(); + fHBlindPixIdTheta->Write(); + + fHSigmaTheta->Write(); + fHSigmaThetaON->Write(); + fHSigmaThetaOFF->Write(); + + fHSigmaPixTheta->Write(); + fHSigmaPixThetaON->Write(); + fHSigmaPixThetaOFF->Write(); + + fHDiffPixTheta->Write(); + fHDiffPixThetaON->Write(); + fHDiffPixThetaOFF->Write(); + + fHgON->Write(); + fHgOFF->Write(); + + *fLog << all << "WriteTargetDist; target histograms written onto file " + << namefileout << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set type of data to be padded +// +// this is not necessary if the type of the data can be recognized +// directly from the input +// +// +void MPadONOFF::SetDataType(const char* type) +{ + fType = type; + *fLog << all << "SetDataType(); type of data to be padded : " + << fType << endl; + + return; +} + + +// -------------------------------------------------------------------------- +// +// Set the option for the padding +// +// There are 2 options for the padding : +// +// 1) fPadFlag = 1 : +// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta +// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using +// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 +// (fHDiffPixTheta). +// +// This is the preferred option as it takes into account the +// correlations between the Sigma of a pixel and Sigmabar. +// +// NOT YET TESTED : +// 2) fPadFlag = 2 : +// Generate a pedestal sigma for each pixel using the 3D-histogram +// Theta, pixel no., Sigma (fHSigmaPixTheta). +// +void MPadONOFF::SetPadFlag(Int_t padflag) +{ + fPadFlag = padflag; + *fLog << all << "SetPadFlag(); choose option " << fPadFlag << endl; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers and prepare the histograms +// +Int_t MPadONOFF::PreProcess(MParList *pList) +{ + if ( !fHSigmaTheta || !fHSigmaThetaON || !fHSigmaThetaOFF || + !fHSigmaPixTheta || !fHSigmaPixThetaON || !fHSigmaPixThetaOFF || + !fHDiffPixTheta || !fHDiffPixThetaON || !fHDiffPixThetaOFF || + !fHBlindPixIdTheta || !fHBlindPixNTheta || + !fHgON || !fHgOFF) + { + *fLog << err << "At least one of the histograms needed for the padding is not defined ... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << err << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << err << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar"); + if (!fSigmabar) + { + *fLog << err << "MSigmabar not found... aborting." << endl; + return kFALSE; + } + + fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); + if (!fBlinds) + { + *fLog << err << "MBlindPixels not found... aborting." << endl; + return kFALSE; + } + + if (fType !="ON" && fType !="OFF" && fType !="MC") + { + *fLog << err << "Type of data to be padded not defined... aborting." << endl; + return kFALSE; + } + + + //-------------------------------------------------------------------- + // histograms for checking the padding + // + // plot pedestal sigmas + fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", + 100, 0.0, 5.0, 100, 0.0, 5.0); + fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); + fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); + + // plot no.of photons (before vs. after padding) + fHPhotons = new TH2D("Photons","Photons: after vs.before padding", + 100, -10.0, 90.0, 100, -10, 90); + fHPhotons->SetXTitle("no.of photons before padding"); + fHPhotons->SetYTitle("no.of photons after padding"); + + // plot of added NSB + fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0); + fHNSB->SetXTitle("no.of added NSB photons"); + fHNSB->SetYTitle("no.of pixels"); + + + //-------------------------------------------------------------------- + + fIter = 20; + + memset(fErrors, 0, sizeof(fErrors)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Do the Padding +// idealy the events to be padded should have been generated without NSB +// +Int_t MPadONOFF::Process() +{ + *fLog << all << "Entry Process();" << endl; + + // this is the conversion factor from the number of photons + // to the number of photo electrons + // later to be taken from MCalibration + // fPEpePhoton = PW * LG * QE * 1D + Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; + + //------------------------------------------------ + // add pixels to MCerPhotEvt which are not yet in; + // set their number of photons equal to zero + + UInt_t imaxnumpix = fCam->GetNumPixels(); + + for (UInt_t i=0; iGetNumPixels(); + for (UInt_t j=0; jAddPixel(i, 0.0, (*fPed)[i].GetPedestalRms()); + } + + + + //----------------------------------------- + Int_t rc=0; + Int_t rw=0; + + const UInt_t npix = fEvt->GetNumPixels(); + + //fSigmabar->Calc(*fCam, *fPed, *fEvt); + // *fLog << all << "before padding : " << endl; + //fSigmabar->Print(""); + + + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + // to simulate the situation that before the padding the NSB and + // electronic noise are zero : set Sigma = 0 for all pixels + //for (UInt_t i=0; ioperator[](i); + // Int_t j = pix.GetPixId(); + + // MPedestalPix &ppix = fPed->operator[](j); + // ppix.SetMeanRms(0.0); + //} + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + + //------------------------------------------- + // Calculate average sigma of the event + // + Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); + *fLog << all << "before padding : " << endl; + fSigmabar->Print(""); + Double_t sigbarold = sigbarold_phot * fPEperPhoton; + + + Double_t sigbarold2 = sigbarold*sigbarold; + + + // for MC data : expect sigmabar to be zero before the padding + if (fType == "MC") + { + if (sigbarold_phot > 0) + { + *fLog << warn << "Process(); sigmabar of event to be padded is > 0 : " + << sigbarold_phot << ". Stop event loop " << endl; + // input data should have sigmabar = 0; stop event loop + + rc = 1; + fErrors[rc]++; + return kCONTINUE; + } + } + + const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + *fLog << all << "theta = " << theta << endl; + + Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta); + if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() ) + { + *fLog << warn << "Process(); binNumber out of range : theta, binTheta = " + << theta << ", " << binTheta << "; Skip event " << endl; + // event cannot be padded; skip event + + rc = 2; + fErrors[rc]++; + return kCONTINUE; + } + + //------------------------------------------- + // get number of events in this theta bin + // and number of events in this sigbarold_phot bin + + Double_t nTheta; + Double_t nSigma; + if (fType == "ON") + { + TH1D *hn; + + hn = fHSigmaThetaON->ProjectionY("", binTheta, binTheta, ""); + nTheta = hn->Integral(); + Int_t binSigma = hn->FindBin(sigbarold_phot); + nSigma = hn->GetBinContent(binSigma); + + *fLog << all + << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " + << theta << ", " << sigbarold_phot << ", " << binTheta << ", " + << binSigma << ", " << nTheta << ", " << nSigma << endl; + + delete hn; + } + + else if (fType == "OFF") + { + TH1D *hn; + + hn = fHSigmaThetaOFF->ProjectionY("", binTheta, binTheta, ""); + nTheta = hn->Integral(); + Int_t binSigma = hn->FindBin(sigbarold_phot); + nSigma = hn->GetBinContent(binSigma); + + *fLog << all + << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " + << theta << ", " << sigbarold_phot << ", " << binTheta << ", " + << binSigma << ", " << nTheta << ", " << nSigma << endl; + + delete hn; + } + + else + { + nTheta = 0.0; + nSigma = 0.0; + } + + //------------------------------------------- + // for the current theta, + // generate blind pixels according to the histograms + // fHBlindPixNTheta and fHBlindPixIDTheta + // do this only for MC data + + + if (fType == "MC") + { + + // numBlind is the number of blind pixels in this event + TH1D *nblind; + UInt_t numBlind; + + nblind = fHBlindPixNTheta->ProjectionY("", binTheta, binTheta, ""); + if ( nblind->Integral() == 0.0 ) + { + *fLog << warn << "Process(); projection for Theta bin " + << binTheta << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete nblind; + + rc = 7; + fErrors[rc]++; + return kCONTINUE; + } + else + { + numBlind = (Int_t) (nblind->GetRandom()+0.5); + } + delete nblind; + + //warn Code commented out due no compilation errors on Alpha OSF (tgb) + + // throw the Id of numBlind different pixels in this event + TH1D *hblind; + UInt_t idBlind; + UInt_t listId[npix]; + UInt_t nlist = 0; + Bool_t equal; + + hblind = fHBlindPixIdTheta->ProjectionY("", binTheta, binTheta, ""); + if ( hblind->Integral() > 0.0 ) + for (UInt_t i=0; iGetRandom()+0.5); + equal = kFALSE; + for (UInt_t j=0; jSetPixelBlind(idBlind); + *fLog << all << "idBlind = " << idBlind << endl; + } + fBlinds->SetReadyToSave(); + + delete hblind; + + } + + //****************************************************************** + // has the event to be padded ? + // if yes : to which sigmabar should it be padded ? + // + + Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot); + *fLog << all << "binSig, sigbarold_phot = " << binSig << ", " + << sigbarold_phot << endl; + + Double_t prob; + TH1D *hpad = NULL; + if (fType == "ON") + { + hpad = fHgON->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, ""); + + //Int_t nb = hpad->GetNbinsX(); + //for (Int_t i=1; i<=nb; i++) + //{ + // if (hpad->GetBinContent(i) != 0.0) + // *fLog << all << "i, fHgON = " << i << ", " + // << hpad->GetBinContent(i) << endl; + //} + + if (nSigma != 0.0) + prob = hpad->Integral() * nTheta/nSigma; + else + prob = 0.0; + + *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma + << ", " << prob << endl; + } + else if (fType == "OFF") + { + hpad = fHgOFF->ProjectionZ("zOFF", binTheta, binTheta, binSig, binSig, ""); + if (nSigma != 0.0) + prob = hpad->Integral() * nTheta/nSigma; + else + prob = 0.0; + } + else + prob = 1.0; + + if ( fType != "MC" && + (prob <= 0.0 || gRandom->Uniform() > prob) ) + { + delete hpad; + // event should not be padded + *fLog << all << "event is not padded" << endl; + + rc = 8; + fErrors[rc]++; + return kTRUE; + } + // event should be padded + *fLog << all << "event is padded" << endl; + + + //------------------------------------------- + // for the current theta, generate a sigmabar + // + // for ON/OFF data according to the matrix fHgON/OFF + // for MC data according to the histogram fHSigmaTheta + // + Double_t sigmabar_phot = 0; + Double_t sigmabar = 0; + if (fType == "ON" || fType == "OFF") + { + //Int_t nbinsx = hpad->GetNbinsX(); + //for (Int_t i=1; i<=nbinsx; i++) + //{ + // if (hpad->GetBinContent(i) != 0.0) + // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i) + // << endl; + //} + + sigmabar_phot = hpad->GetRandom(); + sigmabar = sigmabar_phot * fPEperPhoton; + + *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl; + + delete hpad; + } + + else if (fType == "MC") + { + Int_t bincheck = fHSigmaTheta->GetXaxis()->FindBin(theta); + + if (binTheta != bincheck) + { + cout << "The binnings of the 2 histograms are not identical; aborting" + << endl; + return kERROR; + } + + TH1D *hsigma; + + hsigma = fHSigmaTheta->ProjectionY("", binTheta, binTheta, ""); + if ( hsigma->Integral() == 0.0 ) + { + *fLog << warn << "Process(); projection for Theta bin " + << binTheta << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete hsigma; + + rc = 3; + fErrors[rc]++; + return kCONTINUE; + } + else + { + sigmabar_phot = hsigma->GetRandom(); + sigmabar = sigmabar_phot * fPEperPhoton; + + *fLog << all << "Theta, bin number = " << theta << ", " << binTheta + << ", sigmabar_phot = " << sigmabar_phot << endl; + } + delete hsigma; + } + + const Double_t sigmabar2 = sigmabar*sigmabar; + + //------------------------------------------- + + *fLog << all << "Process(); sigbarold, sigmabar = " + << sigbarold << ", "<< sigmabar << endl; + + // Skip event if target sigmabar is <= sigbarold + if (sigmabar <= sigbarold) + { + *fLog << all << "Process(); target sigmabar is less than sigbarold : " + << sigmabar << ", " << sigbarold << ", Skip event" << endl; + + rc = 4; + fErrors[rc]++; + return kCONTINUE; + } + + + //------------------------------------------- + // + // Calculate average number of NSB photo electrons to be added (lambdabar) + // from the value of sigmabar, + // - making assumptions about the average electronic noise (elNoise2) and + // - using a fixed value (F2excess) for the excess noise factor + + Double_t elNoise2; // [photons] + Double_t F2excess = 1.3; + Double_t lambdabar; // [photons] + + + + Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta); + if (binTheta != bincheck) + { + cout << "The binnings of the 2 histograms are not identical; aborting" + << endl; + return kERROR; + } + + // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin. + // The average electronic noise (to be added) has to be well above this RMS, + // otherwise the electronic noise of an individual pixel (elNoise2Pix) + // may become negative + + TH1D *hnoise; + if (fType == "MC") + hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else if (fType == "ON") + hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else if (fType == "OFF") + hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else + { + *fLog << err << "Process; illegal data type... aborting" << endl; + return kERROR; + } + + Double_t RMS_phot = hnoise->GetRMS(1); + Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton; + delete hnoise; + + elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); + *fLog << all << "elNoise2 = " << elNoise2 << endl; + + lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; // [photo electrons] + + // This value of lambdabar is the same for all pixels; + // note that lambdabar is normalized to the area of pixel 0 + + //---------- start loop over pixels --------------------------------- + // do the padding for each pixel + // + // pad only pixels - which are used (before image cleaning) + // + Double_t sig_phot = 0; + Double_t sig = 0; + + Double_t sigma2 = 0; + + Double_t diff_phot = 0; + Double_t diff = 0; + + Double_t addSig2_phot= 0; + Double_t addSig2 = 0; + + Double_t elNoise2Pix = 0; + + + for (UInt_t i=0; iGetPixRatio(j); + + MPedestalPix &ppix = (*fPed)[j]; + Double_t oldsigma_phot = ppix.GetPedestalRms(); + Double_t oldsigma = oldsigma_phot * fPEperPhoton; + Double_t oldsigma2 = oldsigma*oldsigma; + + //--------------------------------- + // throw the Sigma for this pixel + // + Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j ); + + Int_t count; + Bool_t ok; + + TH1D *hdiff; + TH1D *hsig; + + switch (fPadFlag) + { + case 1 : + // throw the Sigma for this pixel from the distribution fHDiffPixTheta + + if (fType == "MC") + hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + else if (fType == "ON") + hdiff = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + else + hdiff = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + + if ( hdiff->Integral() == 0 ) + { + *fLog << warn << "Process(); projection for Theta bin " + << binTheta << " and pixel bin " << binPixel + << " has no entries; aborting " << endl; + delete hdiff; + + rc = 5; + fErrors[rc]++; + return kCONTINUE; + } + + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + diff = diff_phot * fPEperPhoton * fPEperPhoton; + + // the following condition ensures that elNoise2Pix > 0.0 + if ( (diff + sigmabar2 - oldsigma2/ratioArea + - lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + if (!ok) + { + *fLog << all << "theta, j, count, sigmabar, diff = " << theta + << ", " + << j << ", " << count << ", " << sigmabar << ", " + << diff << endl; + diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2; + + rw = 1; + fWarnings[rw]++; + } + else + { + rw = 0; + fWarnings[rw]++; + } + + delete hdiff; + sigma2 = diff + sigmabar2; + break; + + case 2 : + // throw the Sigma for this pixel from the distribution fHSigmaPixTheta + + if (fType == "MC") + hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + else if (fType == "ON") + hsig = fHSigmaPixThetaOFF->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + else + hsig = fHSigmaPixThetaON->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + + if ( hsig->Integral() == 0 ) + { + *fLog << warn << "Process(); projection for Theta bin " + << binTheta << " and pixel bin " << binPixel + << " has no entries; aborting " << endl; + delete hsig; + + rc = 6; + fErrors[rc]++; + return kCONTINUE; + } + + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + sig = sig_phot * fPEperPhoton; + + sigma2 = sig*sig/ratioArea; + // the following condition ensures that elNoise2Pix > 0.0 + + if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + if (!ok) + { + *fLog << all << "theta, j, count, sigmabar, sig = " << theta + << ", " + << j << ", " << count << ", " << sigmabar << ", " + << sig << endl; + sigma2 = lambdabar*F2excess + oldsigma2/ratioArea; + + rw = 1; + fWarnings[rw]++; + } + else + { + rw = 0; + fWarnings[rw]++; + } + + delete hsig; + break; + } + + //--------------------------------- + // get the additional sigma^2 for this pixel (due to the padding) + + addSig2 = sigma2*ratioArea - oldsigma2; + addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); + + //--------------------------------- + // get the additional electronic noise for this pixel + + elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea; + + + //--------------------------------- + // throw actual number of additional NSB photo electrons (NSB) + // and its RMS (sigmaNSB) + + Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea); + Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix; + Double_t sigmaNSB0; + + if (arg >= 0.0) + { + sigmaNSB0 = sqrt( arg ); + } + else + { + sigmaNSB0 = 0.0000001; + if (arg < -1.e-10) + { + *fLog << err << "Process(); argument of sqrt < 0 : " + << arg << endl; + } + } + + + //--------------------------------- + // smear NSB0 according to sigmaNSB0 + // and subtract lambdabar because of AC coupling + + Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea; + Double_t NSB_phot = NSB / fPEperPhoton; + + //--------------------------------- + + // add additional NSB to the number of photons + Double_t oldphotons_phot = pix.GetNumPhotons(); + Double_t oldphotons = oldphotons_phot * fPEperPhoton; + Double_t newphotons = oldphotons + NSB; + Double_t newphotons_phot = newphotons / fPEperPhoton;; + pix.SetNumPhotons( newphotons_phot ); + + + fHNSB->Fill( NSB_phot/ratioArea ); + fHPhotons->Fill( oldphotons_phot/ratioArea, + newphotons_phot/ratioArea ); + + + // error: add sigma of padded noise quadratically + Double_t olderror_phot = pix.GetErrorPhot(); + Double_t newerror_phot = + sqrt( olderror_phot*olderror_phot + addSig2_phot ); + pix.SetErrorPhot( newerror_phot ); + + + Double_t newsigma = sqrt( oldsigma2 + addSig2 ); + Double_t newsigma_phot = newsigma / fPEperPhoton;; + ppix.SetPedestalRms( newsigma_phot ); + + fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot ); + } + //---------- end of loop over pixels --------------------------------- + + // Calculate sigmabar again and crosscheck + + fSigmabar->Calc(*fCam, *fPed, *fEvt); + *fLog << all << "after padding : " << endl; + fSigmabar->Print(""); + + rc = 0; + fErrors[rc]++; + return kTRUE; + //****************************************************************** +} + +// -------------------------------------------------------------------------- +// +// +Int_t MPadONOFF::PostProcess() +{ + if (GetNumExecutions() != 0) + { + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + + int temp; + if (fWarnings[0] != 0.0) + temp = (int)(fWarnings[1]*100/fWarnings[0]); + else + temp = 100; + + *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) + << temp + << "%) Warning : iteration to find acceptable sigma failed" + << ", fIter = " << fIter << endl; + + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) + << (int)(fErrors[1]*100/GetNumExecutions()) + << "%) Evts skipped due to: Sigmabar_old > 0" << endl; + + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) + << (int)(fErrors[2]*100/GetNumExecutions()) + << "%) Evts skipped due to: Zenith angle out of range" << endl; + + *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) + << (int)(fErrors[3]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigmabar" << endl; + + *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) + << (int)(fErrors[4]*100/GetNumExecutions()) + << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; + + *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3) + << (int)(fErrors[5]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl; + + *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) + << (int)(fErrors[6]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigma" << endl; + + *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) + << (int)(fErrors[7]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Blind pixels" << endl; + + *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3) + << (int)(fErrors[8]*100/GetNumExecutions()) + << "%) Evts didn't have to be padded" << endl; + + *fLog << " " << fErrors[0] << " (" + << (int)(fErrors[0]*100/GetNumExecutions()) + << "%) Evts were successfully padded!" << endl; + *fLog << endl; + + } + + //--------------------------------------------------------------- + TCanvas &c = *(MH::MakeDefCanvas("PadONOFF", "", 900, 1200)); + c.Divide(3, 4); + gROOT->SetSelectedPad(NULL); + + c.cd(1); + fHNSB->SetDirectory(NULL); + fHNSB->DrawCopy(); + fHNSB->SetBit(kCanDelete); + + c.cd(2); + fHSigmaPedestal->SetDirectory(NULL); + fHSigmaPedestal->DrawCopy(); + fHSigmaPedestal->SetBit(kCanDelete); + + c.cd(3); + fHPhotons->SetDirectory(NULL); + fHPhotons->DrawCopy(); + fHPhotons->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(4); + fHSigmaTheta->SetDirectory(NULL); + fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta"); + fHSigmaTheta->DrawCopy(); + fHSigmaTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + + + c.cd(7); + fHBlindPixNTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetTitle("(Target) 2D : no.of blind pixels, \\Theta"); + fHBlindPixNTheta->DrawCopy(); + fHBlindPixNTheta->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(10); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetTitle("(Target) 2D : blind pixel Id, \\Theta"); + fHBlindPixIdTheta->DrawCopy(); + fHBlindPixIdTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2 + + c.cd(5); + TH2D *l1; + l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx"); + l1->SetDirectory(NULL); + l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); + l1->SetXTitle("\\Theta [\\circ]"); + l1->SetYTitle("Sigma^2-Sigmabar^2"); + + l1->DrawCopy("box"); + l1->SetBit(kCanDelete);; + + c.cd(8); + TH2D *l2; + l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy"); + l2->SetDirectory(NULL); + l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); + l2->SetXTitle("pixel"); + l2->SetYTitle("Sigma^2-Sigmabar^2"); + + l2->DrawCopy("box"); + l2->SetBit(kCanDelete);; + + //-------------------------------------------------------------------- + // draw the 3D histogram (target): Theta, pixel, Sigma + + c.cd(6); + TH2D *k1; + k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx"); + k1->SetDirectory(NULL); + k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)"); + k1->SetXTitle("\\Theta [\\circ]"); + k1->SetYTitle("Sigma"); + + k1->DrawCopy("box"); + k1->SetBit(kCanDelete); + + c.cd(9); + TH2D *k2; + k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy"); + k2->SetDirectory(NULL); + k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)"); + k2->SetXTitle("pixel"); + k2->SetYTitle("Sigma"); + + k2->DrawCopy("box"); + k2->SetBit(kCanDelete);; + + + //-------------------------------------------------------------------- + + c.cd(11); + TH2D *m1; + m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); + m1->SetDirectory(NULL); + m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); + m1->SetXTitle("Sigmabar-old"); + m1->SetYTitle("Sigmabar-new"); + + m1->DrawCopy("box"); + m1->SetBit(kCanDelete);; + + c.cd(12); + TH2D *m2; + m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); + m2->SetDirectory(NULL); + m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); + m2->SetXTitle("Sigmabar-old"); + m2->SetYTitle("Sigmabar-new"); + + m2->DrawCopy("box"); + m2->SetBit(kCanDelete);; + + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MPadONOFF.h (revision 2663) @@ -0,0 +1,105 @@ +#ifndef MARS_MPadONOFF +#define MARS_MPadONOFF + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TH2D; +class TH3D; + +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MMcEvt; +class MSigmabar; +class MParList; +class MBlindPixels; +class MRead; +class MFilterList; + + +class MPadONOFF : public MTask +{ +private: + MGeomCam *fCam; + MCerPhotEvt *fEvt; + MSigmabar *fSigmabar; + MMcEvt *fMcEvt; + MPedestalCam *fPed; + MBlindPixels *fBlinds; + + TString fType; // type of data to be padded + TFile *fInfile; // input file containing padding histograms + + Int_t fPadFlag; + Int_t fIter; + + Int_t fErrors[9]; + Int_t fWarnings[2]; + + + // plots used for the padding + // for all plots it is assumed that the pedestal RMS is given in units of "number of photons" + TH2D *fHBlindPixIdTheta; // 2D-histogram (blind pixel Id vs. Theta) + TH2D *fHBlindPixNTheta; // 2D-histogram (no.of blind pixels vs. Theta) + + TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) + TH2D *fHSigmaThetaON; // 2D-histogram (sigmabar vs. Theta) + TH2D *fHSigmaThetaOFF; // 2D-histogram (sigmabar vs. Theta) + + TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHSigmaPixThetaON; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHSigmaPixThetaOFF; // 3D-histogram (Theta, pixel, sigma) + + TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + TH3D *fHDiffPixThetaON; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + TH3D *fHDiffPixThetaOFF; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + + TH3D *fHgON; // matrix (Theta, sigbarold, sigbarnew) for ON data + TH3D *fHgOFF; // matrix (Theta, sigbarold, sigbarnew) for OFF data + + // plots for checking the padding + TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after + // versus before padding + TH2D *fHPhotons; // 2D-histogram : no.of photons after + // versus before padding + TH1D *fHNSB; // 1D-histogram : additional NSB + + +public: + MPadONOFF(const char *name=NULL, const char *title=NULL); + ~MPadONOFF(); + + Bool_t MergeHistograms(TH2D *sigthon, TH2D *sigthoff, + TH3D *sigpixthon, TH3D *sigpixthoff, + TH3D *diffpixthon, TH3D *diffpixthoff, + TH2D *blindidthon, TH2D *blindidthoff, + TH2D *blindnthon, TH2D *blindnthoff); + + Bool_t ReadTargetDist(const char *filein); + Bool_t WriteTargetDist(const char *fileout); + + void SetDataType(const char *type); // type of data to be padded + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + void SetPadFlag(Int_t padflag); + + ClassDef(MPadONOFF, 0) // task for the ON-OFF padding +}; + +#endif + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc (revision 2663) @@ -0,0 +1,114 @@ +/* ======================================================================== *\ +! +! * +! * 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): Sabrina Stark 12/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MPedestalWorkaround // +// // +// This class contains the function to store the pedestal value and the // +// RMS of the pedestal in units of photons. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MPedestalWorkaround.h" + +#include +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MGeomCam.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MPedPhotCam.h" + +using namespace std; + +ClassImp(MPedestalWorkaround); + +MPedestalWorkaround::MPedestalWorkaround(const char *name, const char *title) +{ + fName = name ? name : "MPedestalWorkaround"; + fTitle = title ? title : "Storage of pedestal values and RMS in units of photons"; + +} + +// ------------------------------------------------------------------------ +// + +Int_t MPedestalWorkaround::PreProcess(MParList *pList) +{ + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << err << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + fPedPhot = (MPedPhotCam*)pList->FindObject("MPedPhotCam"); + if (!fPedPhot) + { + *fLog << err << "MPedPhotCam not found... aborting." << endl; + return kFALSE; + } + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// ------------------------------------------------------------------------ +// +Int_t MPedestalWorkaround::Process() +{ + + UInt_t imaxnumpix = fCam->GetNumPixels(); + + for (UInt_t i=0; iGetPixelContent( val, i, *fCam, type); + valout = (*fPed)[i].GetPedestal(); + (*fPed)[i].SetPedestal(val); + *fLog << "i, val, valout : " << i <<", "<< val<<", " << valout << endl; + type = 1; + fPedPhot->GetPixelContent( val, i, *fCam, type); + valout = (*fPed)[i].GetPedestalRms(); + (*fPed)[i].SetPedestalRms(val); + *fLog << "RMS : i, val, valout : " << i <<", "<< val<<", " << valout << endl; + } + + + return kTRUE; +} + +Int_t MPedestalWorkaround::PostProcess() +{ + return kTRUE; +} + Index: trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h (revision 2663) @@ -0,0 +1,44 @@ +#ifndef MARS_MPedestalWorkaround +#define MARS_MPedestalWorkaround + +///////////////////////////////////////////////////////////////////////////// +// // +// MPedestalWorkaround // +// // +// Temporary copy of pedestal values and RMS into MPedestalCam from // +// MPedPhotCam in order to have the pedestal RMS in units of number of // +// photons // +// // +///////////////////////////////////////////////////////////////////////////// +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MPedestalCam; +class MPedPhotCam; +class MGeomCam; +class MParList; + + +class MPedestalWorkaround : public MTask +{ +private: + + MPedestalCam *fPed; // + MPedPhotCam *fPedPhot; // + MGeomCam *fCam; // + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + MPedestalWorkaround(const char *name=NULL, const char *title=NULL); + + ClassDef(MPedestalWorkaround, 0) // Task to copy pedestal RMS +}; + +#endif + + + Index: trunk/MagicSoft/Mars/manalysis/MSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercuts.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MSupercuts.cc (revision 2663) @@ -0,0 +1,388 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek, 08/2003 +! Author(s): Thomas Bretz, 08/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSupercuts // +// // +// this is the container for the parameters of the supercuts // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MSupercuts.h" + +#include +#include + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSupercuts); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// constructor +// +MSupercuts::MSupercuts(const char *name, const char *title) + : fParameters(104), fStepsizes(104), + fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+8), + fWidthUp(fParameters.GetArray()+16), fWidthLo(fParameters.GetArray()+24), + fDistUp(fParameters.GetArray()+32), fDistLo(fParameters.GetArray()+40), + fAsymUp(fParameters.GetArray()+48), fAsymLo(fParameters.GetArray()+56), + fConcUp(fParameters.GetArray()+64), fConcLo(fParameters.GetArray()+72), + fLeakage1Up(fParameters.GetArray()+80), fLeakage1Lo(fParameters.GetArray()+88), + fAlphaUp(fParameters.GetArray()+96) +{ + fName = name ? name : "MSupercuts"; + fTitle = title ? title : "Container for the supercut parameters"; + + // set supercut parameters to their default values + InitParameters(); +} + + +// -------------------------------------------------------------------------- +// +// set default values for the supercut parameters +// +void MSupercuts::InitParameters() +{ + //--------------------------------------------------- + // these are Daniel's original values for Mkn 421 + + fLengthUp[0] = 0.315585; + fLengthUp[1] = 0.001455; + fLengthUp[2] = 0.203198; + fLengthUp[3] = 0.005532; + fLengthUp[4] = -0.001670; + fLengthUp[5] = -0.020362; + fLengthUp[6] = 0.007388; + fLengthUp[7] = -0.013463; + + fLengthLo[0] = 0.151530; + fLengthLo[1] = 0.028323; + fLengthLo[2] = 0.510707; + fLengthLo[3] = 0.053089; + fLengthLo[4] = 0.013708; + fLengthLo[5] = 2.357993; + fLengthLo[6] = 0.000080; + fLengthLo[7] = -0.007157; + + fWidthUp[0] = 0.145412; + fWidthUp[1] = -0.001771; + fWidthUp[2] = 0.054462; + fWidthUp[3] = 0.022280; + fWidthUp[4] = -0.009893; + fWidthUp[5] = 0.056353; + fWidthUp[6] = 0.020711; + fWidthUp[7] = -0.016703; + + fWidthLo[0] = 0.089187; + fWidthLo[1] = -0.006430; + fWidthLo[2] = 0.074442; + fWidthLo[3] = 0.003738; + fWidthLo[4] = -0.004256; + fWidthLo[5] = -0.014101; + fWidthLo[6] = 0.006126; + fWidthLo[7] = -0.002849; + + fDistUp[0] = 1.787943; + fDistUp[1] = 0; + fDistUp[2] = 2.942310; + fDistUp[3] = 0.199815; + fDistUp[4] = 0; + fDistUp[5] = 0.249909; + fDistUp[6] = 0.189697; + fDistUp[7] = 0; + + fDistLo[0] = 0.589406; + fDistLo[1] = 0; + fDistLo[2] = -0.083964; + fDistLo[3] = -0.007975; + fDistLo[4] = 0; + fDistLo[5] = 0.045374; + fDistLo[6] = -0.001750; + fDistLo[7] = 0; + + + // dummy values + + fAsymUp[0] = 1.e10; + fAsymUp[1] = 0.0; + fAsymUp[2] = 0.0; + fAsymUp[3] = 0.0; + fAsymUp[4] = 0.0; + fAsymUp[5] = 0.0; + fAsymUp[6] = 0.0; + fAsymUp[7] = 0.0; + + fAsymLo[0] = -1.e10; + fAsymLo[1] = 0.0; + fAsymLo[2] = 0.0; + fAsymLo[3] = 0.0; + fAsymLo[4] = 0.0; + fAsymLo[5] = 0.0; + fAsymLo[6] = 0.0; + fAsymLo[7] = 0.0; + + fConcUp[0] = 1.e10; + fConcUp[1] = 0.0; + fConcUp[2] = 0.0; + fConcUp[3] = 0.0; + fConcUp[4] = 0.0; + fConcUp[5] = 0.0; + fConcUp[6] = 0.0; + fConcUp[7] = 0.0; + + fConcLo[0] = -1.e10; + fConcLo[1] = 0.0; + fConcLo[2] = 0.0; + fConcLo[3] = 0.0; + fConcLo[4] = 0.0; + fConcLo[5] = 0.0; + fConcLo[6] = 0.0; + fConcLo[7] = 0.0; + + fLeakage1Up[0] = 1.e10; + fLeakage1Up[1] = 0.0; + fLeakage1Up[2] = 0.0; + fLeakage1Up[3] = 0.0; + fLeakage1Up[4] = 0.0; + fLeakage1Up[5] = 0.0; + fLeakage1Up[6] = 0.0; + fLeakage1Up[7] = 0.0; + + fLeakage1Lo[0] = -1.e10; + fLeakage1Lo[1] = 0.0; + fLeakage1Lo[2] = 0.0; + fLeakage1Lo[3] = 0.0; + fLeakage1Lo[4] = 0.0; + fLeakage1Lo[5] = 0.0; + fLeakage1Lo[6] = 0.0; + fLeakage1Lo[7] = 0.0; + + fAlphaUp[0] = 13.123440; + fAlphaUp[1] = 0; + fAlphaUp[2] = 0; + fAlphaUp[3] = 0; + fAlphaUp[4] = 0; + fAlphaUp[5] = 0; + fAlphaUp[6] = 0; + fAlphaUp[7] = 0; + + //--------------------------------------------------- + // fStepsizes + // if == 0.0 the parameter will be fixed in the minimization + // != 0.0 initial step sizes for the parameters + + // LengthUp + fStepsizes[0] = 0.03; + fStepsizes[1] = 0.0002; + fStepsizes[2] = 0.02; + fStepsizes[3] = 0.0006; + fStepsizes[4] = 0.0002; + fStepsizes[5] = 0.002; + fStepsizes[6] = 0.0008; + fStepsizes[7] = 0.002; + + // LengthLo + fStepsizes[8] = 0.02; + fStepsizes[9] = 0.003; + fStepsizes[10] = 0.05; + fStepsizes[11] = 0.006; + fStepsizes[12] = 0.002; + fStepsizes[13] = 0.3; + fStepsizes[14] = 0.0001; + fStepsizes[15] = 0.0008; + + // WidthUp + fStepsizes[16] = 0.02; + fStepsizes[17] = 0.0002; + fStepsizes[18] = 0.006; + fStepsizes[19] = 0.003; + fStepsizes[20] = 0.002; + fStepsizes[21] = 0.006; + fStepsizes[22] = 0.002; + fStepsizes[23] = 0.002; + + // WidthLo + fStepsizes[24] = 0.009; + fStepsizes[25] = 0.0007; + fStepsizes[26] = 0.008; + fStepsizes[27] = 0.0004; + fStepsizes[28] = 0.0005; + fStepsizes[29] = 0.002; + fStepsizes[30] = 0.0007; + fStepsizes[31] = 0.003; + + // DistUp + fStepsizes[32] = 0.2; + fStepsizes[33] = 0.0; + fStepsizes[34] = 0.3; + fStepsizes[35] = 0.02; + fStepsizes[36] = 0.0; + fStepsizes[37] = 0.03; + fStepsizes[38] = 0.02; + fStepsizes[39] = 0.0; + + // DistLo + fStepsizes[40] = 0.06; + fStepsizes[41] = 0.0; + fStepsizes[42] = 0.009; + fStepsizes[43] = 0.0008; + fStepsizes[44] = 0.0; + fStepsizes[45] = 0.005; + fStepsizes[46] = 0.0002; + fStepsizes[47] = 0.0; + + // AsymUp + fStepsizes[48] = 0.0; + fStepsizes[49] = 0.0; + fStepsizes[50] = 0.0; + fStepsizes[51] = 0.0; + fStepsizes[52] = 0.0; + fStepsizes[53] = 0.0; + fStepsizes[54] = 0.0; + fStepsizes[55] = 0.0; + + // AsymLo + fStepsizes[56] = 0.0; + fStepsizes[57] = 0.0; + fStepsizes[58] = 0.0; + fStepsizes[59] = 0.0; + fStepsizes[60] = 0.0; + fStepsizes[61] = 0.0; + fStepsizes[62] = 0.0; + fStepsizes[63] = 0.0; + + // ConcUp + fStepsizes[64] = 0.0; + fStepsizes[65] = 0.0; + fStepsizes[66] = 0.0; + fStepsizes[67] = 0.0; + fStepsizes[68] = 0.0; + fStepsizes[69] = 0.0; + fStepsizes[70] = 0.0; + fStepsizes[71] = 0.0; + + // ConcLo + fStepsizes[72] = 0.0; + fStepsizes[73] = 0.0; + fStepsizes[74] = 0.0; + fStepsizes[75] = 0.0; + fStepsizes[76] = 0.0; + fStepsizes[77] = 0.0; + fStepsizes[78] = 0.0; + fStepsizes[79] = 0.0; + + // Leakage1Up + fStepsizes[80] = 0.0; + fStepsizes[81] = 0.0; + fStepsizes[82] = 0.0; + fStepsizes[83] = 0.0; + fStepsizes[84] = 0.0; + fStepsizes[85] = 0.0; + fStepsizes[86] = 0.0; + fStepsizes[87] = 0.0; + + // Leakage1Lo + fStepsizes[88] = 0.0; + fStepsizes[89] = 0.0; + fStepsizes[90] = 0.0; + fStepsizes[91] = 0.0; + fStepsizes[92] = 0.0; + fStepsizes[93] = 0.0; + fStepsizes[94] = 0.0; + fStepsizes[95] = 0.0; + + // AlphaUp + fStepsizes[96] = 0.0; + fStepsizes[97] = 0.0; + fStepsizes[98] = 0.0; + fStepsizes[99] = 0.0; + fStepsizes[100] = 0.0; + fStepsizes[101] = 0.0; + fStepsizes[102] = 0.0; + fStepsizes[103] = 0.0; +} + + +// -------------------------------------------------------------------------- +// +// Set the parameter values from the array 'd' +// +// +Bool_t MSupercuts::SetParameters(const TArrayD &d) +{ + if (d.GetSize() != fParameters.GetSize()) + { + *fLog << err << "Sizes of d and of fParameters are different : " + << d.GetSize() << ", " << fParameters.GetSize() << endl; + return kFALSE; + } + + fParameters = d; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set the step sizes from the array 'd' +// +// +Bool_t MSupercuts::SetStepsizes(const TArrayD &d) +{ + if (d.GetSize() != fStepsizes.GetSize()) + { + *fLog << err << "Sizes of d and of fStepsizes are different : " + << d.GetSize() << ", " << fStepsizes.GetSize() << endl; + return kFALSE; + } + + fStepsizes = d; + + return kTRUE; +} + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MSupercuts.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercuts.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MSupercuts.h (revision 2663) @@ -0,0 +1,69 @@ +#ifndef MARS_MSupercuts +#define MARS_MSupercuts + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MSupercuts : public MParContainer +{ +private: + TArrayD fParameters; // supercut parameters + TArrayD fStepsizes; // step sizes of supercut parameters + + Double_t *fLengthUp; //! + Double_t *fLengthLo; //! + Double_t *fWidthUp; //! + Double_t *fWidthLo; //! + Double_t *fDistUp; //! + Double_t *fDistLo; //! + Double_t *fAsymUp; //! + Double_t *fAsymLo; //! + + Double_t *fConcUp; //! + Double_t *fConcLo; //! + Double_t *fLeakage1Up; //! + Double_t *fLeakage1Lo; //! + + Double_t *fAlphaUp; //! + + +public: + MSupercuts(const char *name=NULL, const char *title=NULL); + + void InitParameters(); + + Bool_t SetParameters(const TArrayD &d); + Bool_t SetStepsizes(const TArrayD &d); + + const TArrayD &GetParameters() const { return fParameters; } + const TArrayD &GetStepsizes() const { return fStepsizes; } + + const Double_t *GetLengthUp() const { return fLengthUp; } + const Double_t *GetLengthLo() const { return fLengthLo; } + const Double_t *GetWidthUp() const { return fWidthUp; } + const Double_t *GetWidthLo() const { return fWidthLo; } + const Double_t *GetDistUp() const { return fDistUp; } + const Double_t *GetDistLo() const { return fDistLo; } + const Double_t *GetAsymUp() const { return fAsymUp; } + const Double_t *GetAsymLo() const { return fAsymLo; } + + const Double_t *GetConcUp() const { return fConcUp; } + const Double_t *GetConcLo() const { return fConcLo; } + + const Double_t *GetLeakage1Up() const { return fLeakage1Up; } + const Double_t *GetLeakage1Lo() const { return fLeakage1Lo; } + + const Double_t *GetAlphaUp() const { return fAlphaUp; } + + ClassDef(MSupercuts, 1) // A container for the Supercut parameters +}; + +#endif + + + Index: trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.cc (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.cc (revision 2663) @@ -0,0 +1,348 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek, 04/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSupercutsCalc // +// // +// this class calculates the hadronness for the supercuts // +// the parameters of the supercuts are taken // +// from the container MSupercuts // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MSupercutsCalc.h" + +#include +#include + +#include "TFile.h" +#include "TArrayD.h" + +#include "MParList.h" +#include "MHillasExt.h" +#include "MHillasSrc.h" +#include "MNewImagePar.h" +#include "MMcEvt.hxx" +#include "MCerPhotEvt.h" +#include "MGeomCam.h" +#include "MHadronness.h" +#include "MHMatrix.h" +#include "MSupercuts.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSupercutsCalc); + +using namespace std; + + +// -------------------------------------------------------------------------- +// +// constructor +// + +MSupercutsCalc::MSupercutsCalc(const char *hilname, + const char *hilsrcname, + const char *name, const char *title) + : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname), + fHilExtName("MHillasExt"), fNewParName("MNewImagePar"), + fSuperName("MSupercuts") +{ + fName = name ? name : "MSupercutsCalc"; + fTitle = title ? title : "Class to evaluate the Supercuts"; + + fMatrix = NULL; +} + +// -------------------------------------------------------------------------- +// +Int_t MSupercutsCalc::PreProcess(MParList *pList) +{ + MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!cam) + { + *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = cam->GetConvMm2Deg(); + + fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName); + if (!fHadronness) + { + *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl; + return kFALSE; + } + + fSuper = (MSupercuts*)pList->FindObject(fSuperName, "MSupercuts"); + if (!fSuper) + { + *fLog << err << fSuperName << " [MSupercuts] not found... aborting." << endl; + return kFALSE; + } + + if (fMatrix) + return kTRUE; + + //----------------------------------------------------------- + fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); + if (!fHil) + { + *fLog << err << fHilName << " [MHillas] not found... aborting." << endl; + return kFALSE; + } + + fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt"); + if (!fHilExt) + { + *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl; + return kFALSE; + } + + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); + if (!fHilSrc) + { + *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; + return kFALSE; + } + + fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar"); + if (!fNewPar) + { + *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculation of upper and lower limits +// +Double_t MSupercutsCalc::CtsMCut(const Double_t* a, Double_t ls, Double_t ct, + Double_t ls2, Double_t dd2) const +{ + // define cut-function + // + // dNOMLOGSIZE = 4.1 (=log(60.0) + // dNOMCOSZA = 1.0 + // + // a: array of cut parameters + // ls: log(SIZE) - dNOMLOGSIZE + // ls2: ls^2 + // ct: Cos(ZA.) - dNOMCOSZA + // dd2: DIST^2 + const Double_t limit = + a[0] + a[1] * dd2 + a[2] * ct + + ls * (a[3] + a[4] * dd2 + a[5] * ct) + + ls2 * (a[6] + a[7] * dd2); + + //*fLog << "MSupercutsCalc::CtsMCut; *a = " + // << *a << ", " << *(a+1) << ", " << *(a+2) << ", " + // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", " + // << *(a+6) << ", " << *(a+7) << endl; + + //*fLog << "MSupercutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls + // << ", " << ls2 << ", " << ct << ", " << dd2 << ", " + // << limit << endl; + + return limit; +} + +// -------------------------------------------------------------------------- +// +// Returns the mapped value from the Matrix +// +Double_t MSupercutsCalc::GetVal(Int_t i) const +{ + + Double_t val = (*fMatrix)[fMap[i]]; + + + //*fLog << "MSupercutsCalc::GetVal; i, fMatrix, fMap, val = " + // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl; + + + return val; +} + +// -------------------------------------------------------------------------- +// +// 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) MEnergyEstParam::Process +// will take the values from the matrix instead of the containers. +// +void MSupercutsCalc::InitMapping(MHMatrix *mat) +{ + if (fMatrix) + return; + + fMatrix = mat; + + fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta"); + fMap[1] = fMatrix->AddColumn("MHillas.fWidth"); + fMap[2] = fMatrix->AddColumn("MHillas.fLength"); + fMap[3] = fMatrix->AddColumn("MHillas.fSize"); + fMap[4] = fMatrix->AddColumn("MHillas.fMeanX"); + fMap[5] = fMatrix->AddColumn("MHillas.fMeanY"); + fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist"); + fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)"); + fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); + fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc"); + fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1"); +} + +// --------------------------------------------------------------------------- +// +// Evaluate dynamical supercuts +// +// set hadronness to 0.25 if cuts are fullfilled +// 0.75 otherwise +// +Int_t MSupercutsCalc::Process() +{ + const Double_t kNomLogSize = 4.1; + const Double_t kNomCosZA = 1.0; + + const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta(); + const Double_t width0 = fMatrix ? GetVal(1) : fHil->GetWidth(); + const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength(); + const Double_t size = fMatrix ? GetVal(3) : fHil->GetSize(); + const Double_t meanx = fMatrix ? GetVal(4) : fHil->GetMeanX(); + const Double_t meany = fMatrix ? GetVal(5) : fHil->GetMeanY(); + const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist(); + + Double_t help; + if (!fMatrix) + help = TMath::Sign(fHilExt->GetM3Long(), + fHilSrc->GetCosDeltaAlpha()); + const Double_t asym0 = fMatrix ? GetVal(8) : help; + const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc(); + const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1(); + + const Double_t newdist = dist0 * fMm2Deg; + + const Double_t dist2 = meanx*meanx + meany*meany; + const Double_t dist = sqrt(dist2) * fMm2Deg; + const Double_t dd2 = dist*dist; + + + const Double_t dmls = log(size) - kNomLogSize; + const Double_t dmls2 = dmls * dmls; + + const Double_t dmcza = cos(theta) - kNomCosZA; + + const Double_t length = length0 * fMm2Deg; + const Double_t width = width0 * fMm2Deg; + const Double_t asym = asym0 * fMm2Deg; + + /* + *fLog << "newdist, length, width, asym, dist, conc, leakage = " + << newdist << ", " << length << ", " << width << ", " + << asym << ", " << dist << ", " << conc << ", " << leakage + << endl; + + *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = " + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLeakage1Up(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLeakage1Lo(), dmls, dmcza, dmls2, dd2) << ", " + << endl; + */ + + + if ( + //dist < 1.05 && + //newdist < 1.05 && + + newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && + newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && + + length < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) && + length > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) && + + width < CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) && + width > CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) && + + asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) && + asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) && + + dist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && + dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && + + conc < CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) && + conc > CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) && + + leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) && + leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2) ) + + fHadronness->SetHadronness(0.25); + else + fHadronness->SetHadronness(0.75); + + //*fLog << "SChadroness = " << fHadronness->GetHadronness() << endl; + + fHadronness->SetReadyToSave(); + + return kTRUE; +} + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.h (revision 2663) +++ trunk/MagicSoft/Mars/manalysis/MSupercutsCalc.h (revision 2663) @@ -0,0 +1,82 @@ +#ifndef MARS_MSupercutsCalc +#define MARS_MSupercutsCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MParList; +class MHillas; +class MHillasSrc; +class MHillasExt; +class MNewImagePar; +class MMcEvt; +class MCerPhotEvt; +class MGeomCam; +class MHadronness; +class MHMatrix; +class MSupercuts; + +class MSupercutsCalc : public MTask +{ +private: + MHillas *fHil; + MHillasSrc *fHilSrc; + MHillasExt *fHilExt; + MNewImagePar *fNewPar; + MMcEvt *fMcEvt; + MHadronness *fHadronness; //! output container for hadronness + MSupercuts *fSuper; // container for supercut parameters + + TString fHadronnessName; // name of container to store hadronness + TString fHilName; + TString fHilSrcName; + TString fHilExtName; + TString fNewParName; + TString fSuperName; // name of container for supercut parameters + + Double_t fMm2Deg; //! + + Int_t fMap[11]; //! + MHMatrix *fMatrix; //! + + Int_t PreProcess(MParList *pList); + Int_t Process(); + + Double_t GetVal(Int_t i) const; + + Double_t CtsMCut(const Double_t* a, Double_t ls, Double_t ct, + Double_t ls2, Double_t dd2) const; + +public: + MSupercutsCalc(const char *hilname="MHillas", + const char *hilsrcname="MHillasSrc", + const char *name=NULL, const char *title=NULL); + + void SetHadronnessName(const TString name) { fHadronnessName = name; } + TString GetHadronnessName() const { return fHadronnessName; } + + void InitMapping(MHMatrix *mat); + void StopMapping() { InitMapping(NULL); } + + + ClassDef(MSupercutsCalc, 0) // A class to evaluate the Supercuts +}; + +#endif + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2662) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2663) @@ -68,7 +68,11 @@ MCT1SupercutsCalc.cc \ MCT1FindSupercuts.cc \ + MSupercuts.cc \ + MSupercutsCalc.cc \ + MFindSupercuts.cc \ MMinuitInterface.cc \ MFiltercutsCalc.cc \ MCT1PadONOFF.cc \ + MPadONOFF.cc \ MCalibrationCalc.cc \ MCalibrationPix.cc \ @@ -76,4 +80,5 @@ MCalibrationPINDiode.cc \ MCalibrationCam.cc \ + MPedestalWorkaround.cc \ MExtractedSignalCam.cc \ MExtractedSignalPix.cc \ @@ -82,4 +87,5 @@ MArrivalTimeCalc.cc \ MMcCalibrationCalc.cc +>>>>>>> 1.65 SRCS = $(SRCFILES)