Index: trunk/MagicSoft/Mars/manalysis/MEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MEnergyEst.cc (revision 6973) +++ (revision ) @@ -1,57 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Thomas Bretz 1/2002 -! Author(s): Wolfgang Wittek 1/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MEnergyEst // -// // -// Storage Container for the estimated energy // -// // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MEnergyEst.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MEnergyEst); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MEnergyEst::MEnergyEst(const char *name, const char *title) : fEnergy(0) -{ - fName = name ? name : "MEnergyEst"; - fTitle = title ? title : "Storage container for the estimated energy [GeV]"; -} - -void MEnergyEst::Print(Option_t *o) const -{ - *fLog << all << GetDescriptor() << ": Eest = " << (int)fEnergy << "GeV \t Ir = " << (int)(fImpact/100) << "m" << endl; -} - Index: trunk/MagicSoft/Mars/manalysis/MEnergyEst.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MEnergyEst.h (revision 6973) +++ (revision ) @@ -1,28 +1,0 @@ -#ifndef MARS_MEnergyEst -#define MARS_MEnergyEst - -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif - -class MEnergyEst : public MParContainer -{ -private: - Double_t fEnergy; // [GeV] Estimated Energy - Double_t fImpact; // [cm] Estimated Impact - -public: - MEnergyEst(const char *name=NULL, const char *title=NULL); - - void SetEnergy(Double_t e) { fEnergy = e; } - void SetImpact(Double_t i) { fImpact = i; } - Double_t GetEnergy() const { return fEnergy; } - Double_t GetImpact() const { return fImpact; } - - void Print(Option_t *o="") const; - - ClassDef(MEnergyEst, 1) // Storage Container for the estimated Energy -}; - -#endif - Index: trunk/MagicSoft/Mars/manalysis/MFindSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MFindSupercuts.cc (revision 6973) +++ (revision ) @@ -1,1153 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 "MFEventSelector.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"; - - //--------------------------- - // 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; -} - -// -------------------------------------------------------------------------- -// -// 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 (!hreftrain.GetHist().GetEntries()==0) - { - *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"); - - 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 (!hreftest.GetHist().GetEntries()==0) - { - *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"); - - 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 (!href.GetHist().GetEntries()==0) - { - *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(); - - MContinue cont(&selector); - cont.SetName("ContTrainTest"); - - Double_t prob = ( (Double_t) howmanytrain ) - / ( (Double_t)(howmanytrain+howmanytest) ); - MFEventSelector split; - split.SetSelectionRatio(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 6973) +++ (revision ) @@ -1,132 +1,0 @@ -#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); - - 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/MHPedestalCam.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHPedestalCam.cc (revision 6973) +++ (revision ) @@ -1,581 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 12/2000 -! Author(s): Markus Gaug 02/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MHPedestalCam -// -// Fills the extracted pedestals of MExtractedSignalCam into -// the MHCalibrationPix-classes MHPedestalPix for every: -// -// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray -// or MHCalibrationCam::fHiGainArray, respectively, depending if -// MArrivalTimePix::IsLoGainUsed() is set. -// -// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), -// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and -// MHCalibrationCam::fAverageHiGainAreas -// -// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), -// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors -// and MHCalibrationCam::fAverageHiGainSectors -// -// Every pedestal is directly taken from MExtractedSignalCam, filled by the -// appropriate extractor. -// -// The pedestals are filled into a histogram and an array, in order to perform -// a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an -// event-by-event basis and written into the corresponding average pixels. -// -// The histograms are fitted to a Gaussian, mean and sigma with its errors -// and the fit probability are extracted. If none of these values are NaN's and -// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%), -// the fit is declared valid. -// Otherwise, the fit is repeated within ranges of the previous mean -// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit()) -// In case this does not make the fit valid, the histogram means and RMS's are -// taken directly (see MHCalibrationPix::BypassFit()). -// -// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas -// from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup) -// -// The number of summed FADC slices is taken to re-normalize -// the result to a pedestal per pixel with the formulas (see MHPedestalPix::Renorm()): -// - Mean Pedestal / slice = Mean Pedestal / Number slices -// - Mean Pedestal Error / slice = Mean Pedestal Error / Number slices -// - Sigma Pedestal / slice = Sigma Pedestal / Sqrt (Number slices) -// - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices) -// -// The class also fills arrays with the signal vs. event number, creates a fourier -// spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the -// projected fourier components follow an exponential distribution. -// -// This same procedure is performed for the average pixels. -// -// The following results are written into MPedestalCam: -// -// - MCalibrationPix::SetMean() -// - MCalibrationPix::SetMeanErr() -// - MCalibrationPix::SetSigma() -// - MCalibrationPix::SetSigmaErr() -// - MCalibrationPix::SetProb() -// - MCalibrationPix::SetNumPickup() -// -// For all averaged areas, the fitted sigma is multiplied with the square root of -// the number involved pixels in order to be able to compare it to the average of -// sigmas in the camera. -// -///////////////////////////////////////////////////////////////////////////// -#include "MHPedestalCam.h" -#include "MHPedestalPix.h" - -#include "MBadPixelsCam.h" -#include "MBadPixelsPix.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" - -#include "MExtractedSignalCam.h" -#include "MExtractedSignalPix.h" - -#include "MPedestalCam.h" -#include "MPedestalPix.h" - -#include "MGeomCam.h" -#include "MGeomPix.h" - -#include "MCalibrationPedCam.h" - -#include "TH1.h" - -ClassImp(MHPedestalCam); - -using namespace std; -// -------------------------------------------------------------------------- -// -// Default constructor. -// -// Initializes: -// - fExtractHiGainSlices to 0. -// - fExtractLoGainSlices to 0. -// - the event frequency to 1200 Hz. -// - the fRenorm flag to kTRUE -// -MHPedestalCam::MHPedestalCam(const char *name, const char *title) - : fExtractHiGainSlices(0.), fExtractLoGainSlices(0.) -{ - - fName = name ? name : "MHPedestalCam"; - fTitle = title ? title : ""; - - SetPulserFrequency(1200); - SetRenorm(); - -} - - - -// -------------------------------------------------------------------------- -// -// Searches pointer to: -// - MPedestalCam -// - MExtractedSignalCam -// -// Searches or creates: -// - MCalibrationPedCam -// -// Retrieves from MExtractedSignalCam: -// - number of used High Gain FADC slices (MExtractedSignalCam::GetNumUsedHiGainFADCSlices()) -// - number of used Low Gain FADC slices (MExtractedSignalCam::GetNumUsedLoGainFADCSlices()) -// -// Initializes, if empty to MGeomCam::GetNumPixels(): -// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray -// -// Initializes, if empty to MGeomCam::GetNumAreas() for: -// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas -// -// Initializes, if empty to MGeomCam::GetNumSectors() for: -// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors -// -// Calls MHCalibrationCam::InitPedHists() for every entry in: -// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray -// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas -// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors -// -// Sets Titles and Names for the Histograms -// - MHCalibrationCam::fAverageHiGainAreas -// - MHCalibrationCam::fAverageHiGainSectors -// -Bool_t MHPedestalCam::ReInitHists(MParList *pList) -{ - - MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); - if (!signal) - { - gLog << err << "Cannot find MExtractedSignalCam... abort." << endl; - return kFALSE; - } - - - fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); - - if (!fPedestals) - { - gLog << err << "Cannot find MPedestalCam ... abort." << endl; - return kFALSE; - } - - const Int_t npixels = fGeom->GetNumPixels(); - const Int_t nsectors = fGeom->GetNumSectors(); - const Int_t nareas = fGeom->GetNumAreas(); - - fCam = (MCalibrationCam*)pList->FindObject("MCalibrationPedCam"); - if (!fCam) - { - fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationPedCam")); - if (!fCam) - return kFALSE; - - fCam->Init(*fGeom); - } - - Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices(); - Float_t sliceslo = signal->GetNumUsedLoGainFADCSlices(); - - if (sliceshi == 0.) - { - gLog << err << "Number of used signal slices in MExtractedSignalCam is zero ... abort." - << endl; - return kFALSE; - } - - if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices ) - { - gLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam ... abort." - << endl; - return kFALSE; - } - - if (fExtractLoGainSlices != 0. && sliceslo != fExtractLoGainSlices ) - { - gLog << err << "Number of used Low Gain signal slices changed in MExtractedSignalCam ... abort." - << endl; - return kFALSE; - } - - fExtractHiGainSlices = sliceshi; - fExtractLoGainSlices = sliceslo; - - if (fHiGainArray->GetEntries()==0) - { - fHiGainArray->Expand(npixels); - for (Int_t i=0; iGetEntries()==0) - { - fLoGainArray->Expand(npixels); - for (Int_t i=0; iGetEntries()==0) - { - fAverageHiGainAreas->Expand(nareas); - - for (Int_t j=0; jSetTitle("Pedestals average Area Idx "); - - InitPedHists((MHPedestalPix&)GetAverageHiGainArea(j),j,fExtractHiGainSlices); - - } - } - - if (fAverageLoGainAreas->GetEntries()==0) - { - fAverageLoGainAreas->Expand(nareas); - - for (Int_t j=0; jSetTitle("Pedestals average Area Idx "); - - InitPedHists((MHPedestalPix&)GetAverageLoGainArea(j),j,fExtractLoGainSlices); - - } - } - - if (fAverageHiGainSectors->GetEntries()==0) - { - fAverageHiGainSectors->Expand(nsectors); - - for (Int_t j=0; jSetTitle("Pedestals average Sector "); - - InitPedHists((MHPedestalPix&)GetAverageHiGainSector(j),j,fExtractHiGainSlices); - - } - } - - if (fAverageLoGainSectors->GetEntries()==0) - { - fAverageLoGainSectors->Expand(nsectors); - - for (Int_t j=0; jSetTitle("Pedestals average Sector "); - - InitPedHists((MHPedestalPix&)GetAverageLoGainSector(j),j,fExtractLoGainSlices); - - } - } - - return kTRUE; -} - - -// ------------------------------------------------------------- -// -// If MBadPixelsPix::IsBad(): -// - calls MHCalibrationPix::SetExcluded() -// -// Calls: -// - MHGausEvents::InitBins() -// - MHCalibrationPix::ChangeHistId(i) -// - MHCalibrationPix::SetEventFrequency(fPulserFrequency) -// - MHPedestalPix::SetNSlices(nslices) -// -void MHPedestalCam::InitPedHists(MHPedestalPix &hist, const Int_t i, const Float_t nslices) -{ - - hist.InitBins(); - hist.ChangeHistId(i); - hist.SetEventFrequency(fPulserFrequency); - - if (fRenorm) - hist.SetNSlices(nslices); - - hist.SetProbLimit(0.); - - TH1F *h = hist.GetHGausHist(); - h->SetTitle( Form("%s%s", h->GetTitle()," Runs: ")); -} -// ------------------------------------------------------------------------------- -// -// Retrieves pointer to MExtractedSignalCam: -// -// Retrieves from MGeomCam: -// - number of pixels -// - number of pixel areas -// - number of sectors -// -// Fills HiGain or LoGain histograms (MHGausEvents::FillHistAndArray()), respectively -// with the signals MExtractedSignalCam::GetExtractedSignalHiGain() and -// MExtractedSignalCam::GetExtractedSignalLoGain(), respectively. -// -Bool_t MHPedestalCam::FillHists(const MParContainer *par, const Stat_t w) -{ - - MExtractedSignalCam *signal = (MExtractedSignalCam*)par; - if (!signal) - { - gLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl; - return kFALSE; - } - - - - const UInt_t npixels = fGeom->GetNumPixels(); - const UInt_t nareas = fGeom->GetNumAreas(); - const UInt_t nsectors = fGeom->GetNumSectors(); - - TArrayF sumareahi(nareas); - TArrayF sumarealo(nareas); - TArrayF sumsectorhi(nsectors); - TArrayD sumsectorlo(nsectors); - TArrayI numareahi(nareas); - TArrayI numarealo(nareas); - TArrayI numsectorhi(nsectors); - TArrayI numsectorlo(nsectors); - - for (UInt_t i=0; iGetSize() <= idx) - return kFALSE; - - if ((*this)[idx].IsExcluded()) - return kFALSE; - - const Float_t ped = (*fPedestals)[idx].GetPedestal(); - const Float_t rms = (*fPedestals)[idx].GetPedestalRms(); - - const Float_t entsqr = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries()); - - const Float_t pederr = rms/entsqr; - const Float_t rmserr = rms/entsqr/2.; - - const Float_t mean = (*this)[idx].GetMean(); - const Float_t meanerr = (*this)[idx].GetMeanErr(); - const Float_t sigma = (*this)[idx].GetSigma() ; - const Float_t sigmaerr = (*this)[idx].GetSigmaErr(); - - switch (type) - { - case 0: - val = mean; - break; - case 1: - val = meanerr; - break; - case 2: - val = sigma; - break; - case 3: - val = sigmaerr; - break; - case 4: - val = (*this)[idx].GetProb(); - break; - case 5: - val = 2.*(mean-ped)/(ped+mean); - break; - case 6: - val = TMath::Sqrt((pederr*pederr + meanerr*meanerr) * (ped*ped + mean*mean)) - *2./(ped+mean)/(ped+mean); - break; - case 7: - val = 2.*(meanerr-pederr)/(pederr + meanerr); - break; - case 8: - val = 2.*(sigma-rms)/(sigma+rms); - break; - case 9: - val = TMath::Sqrt((rmserr*rmserr + sigmaerr*sigmaerr) * (rms*rms + sigma*sigma)) - *2./(rms+sigma)/(rms+sigma); - break; - case 10: - val = 2.*(sigmaerr - rmserr)/(sigmaerr + rmserr); - break; - case 11: - if (!(*this)[idx].IsGausFitOK()) - val = 1.; - break; - case 12: - if (!(*this)[idx].IsFourierSpectrumOK()) - val = 1.; - break; - default: - return kFALSE; - } - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Calls MHGausEvents::DrawClone() for pixel idx -// -void MHPedestalCam::DrawPixelContent(Int_t idx) const -{ - (*this)[idx].DrawClone(); -} Index: trunk/MagicSoft/Mars/manalysis/MHPedestalCam.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHPedestalCam.h (revision 6973) +++ (revision ) @@ -1,42 +1,0 @@ -#ifndef MARS_MHPedestalCam -#define MARS_MHPedestalCam - -#ifndef MARS_MHCalibrationCam -#include "MHCalibrationCam.h" -#endif - -class MGeomCam; -class MPedestalCam; -class MHPedestalPix; -class MHPedestalCam : public MHCalibrationCam -{ - -private: - - void InitPedHists(MHPedestalPix &hist, const Int_t i, const Float_t nslices); - - Bool_t ReInitHists(MParList *pList); - Bool_t FillHists(const MParContainer *par, const Stat_t w=1); - Bool_t FinalizeHists(); - - Float_t fExtractHiGainSlices; // Number of FADC slices used for high gain signal extraction - Float_t fExtractLoGainSlices; // Number of FADC slices used for low gain signal extraction - - Bool_t fRenorm; - - MPedestalCam *fPedestals; //! Pedestal Cam filled by MPedCalcPedRun - -public: - MHPedestalCam(const char *name=NULL, const char *title=NULL); - ~MHPedestalCam() {} - - Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; - void DrawPixelContent(Int_t idx) const; - - void SetRenorm( const Bool_t b=kTRUE ) { fRenorm = b; } - - ClassDef(MHPedestalCam, 1) // Histogram class for Charge Camera Pedestals -}; - -#endif - Index: trunk/MagicSoft/Mars/manalysis/MHPedestalPix.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHPedestalPix.cc (revision 6973) +++ (revision ) @@ -1,134 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Markus Gaug 02/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MHPedestalPix -// -// Histogram class for pedestal analysis. -// Stores and fits the pedestals taken from MPedestalPix on an event-by-event -// basis. The results are re-normalized to a value per slice with the formulae: -// -// - Mean Pedestal / slice = Mean Pedestal / Number slices -// - Mean Pedestal Error / slice = Mean Pedestal Error / Number slices -// - Sigma Pedestal / slice = Sigma Pedestal / Sqrt (Number slices) -// - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices) -// -// Derives from MHCalibrationPix, fits the pedestals to a Gaussian and performs -// a Fourier analysis. -// -////////////////////////////////////////////////////////////////////////////// -#include "MHPedestalPix.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include - -ClassImp(MHPedestalPix); - -using namespace std; -// -const Int_t MHPedestalPix::fgChargeNbins = 500; -const Axis_t MHPedestalPix::fgChargeFirst = -49.5; -const Axis_t MHPedestalPix::fgChargeLast = 449.5; -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -// Sets: -// - the default number for fNbins (fgChargeNbins) -// - the default number for fFirst (fgChargeFirst) -// - the default number for fLast (fgChargeLast) -// -// - the default name of the fHGausHist ("HPedestalCharge") -// - the default title of the fHGausHist ("Distribution of Summed FADC Pedestal slices Pixel ") -// - the default x-axis title for fHGausHist ("Sum FADC Slices") -// - the default y-axis title for fHGausHist ("Nr. of events") -// - TH1::Sumw2() for fHGausHist -// -// Initializes: -// - fNSlices to 1 -// - fProbLimit to 0.01 -// -MHPedestalPix::MHPedestalPix(const char *name, const char *title) - : fNSlices(1.) -{ - - fName = name ? name : "MHPedestalPix"; - fTitle = title ? title : "Histogrammed Pedestal events"; - - SetNbins( fgChargeNbins ); - SetFirst( fgChargeFirst ); - SetLast( fgChargeLast ); - - SetProbLimit(0.01); - - // Create a large number of bins, later we will rebin - fHGausHist.SetName("HPedestalCharge"); - fHGausHist.SetTitle("Distribution of Summed FADC Pedestal Slices Pixel "); - fHGausHist.SetXTitle("Sum FADC Slices"); - fHGausHist.SetYTitle("Nr. of events"); - fHGausHist.Sumw2(); - -} - -// -------------------------------------------------------------------------- -// -// If mean and sigma have not yet been set, returns. -// -// Renormalizes the pedestal fit results by the following formulae: -// -// - Mean Pedestal / slice = Mean Pedestal / Number slices -// - Mean Pedestal Error / slice = Mean Pedestal Error / Number slices -// - Sigma Pedestal / slice = Sigma Pedestal / Sqrt (Number slices) -// - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices) -// -void MHPedestalPix::Renorm() -{ - - // - // One never knows... - // - if (fNSlices <= 0) - return; - - const Float_t sqslices = TMath::Sqrt(fNSlices); - - SetMean ( GetMean() / fNSlices ); - // - // Mean error goes with PedestalRMS/Sqrt(entries) -> scale with slices - // - SetMeanErr ( GetMeanErr() / fNSlices ); - // - // Sigma goes like PedestalRMS -> scale with sqrt(slices) - // - SetSigma ( GetSigma() / sqslices ); - // - // Sigma error goes like PedestalRMS/2.(entries) -> scale with sqrt(slices) - // - SetSigmaErr ( GetSigmaErr() / sqslices ); - -} - Index: trunk/MagicSoft/Mars/manalysis/MHPedestalPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHPedestalPix.h (revision 6973) +++ (revision ) @@ -1,36 +1,0 @@ -#ifndef MARS_MHPedestalPix -#define MARS_MHPedestalPix - -#ifndef MARS_MHCalibrationPix -#include "MHCalibrationPix.h" -#endif - -class MHPedestalPix : public MHCalibrationPix -{ - -private: - - static const Int_t fgChargeNbins; // Default for fNbins (now set to: 450 ) - static const Axis_t fgChargeFirst; // Default for fFirst (now set to: -0.5 ) - static const Axis_t fgChargeLast; // Default for fLast (now set to: 449.5) - - Float_t fNSlices; // Number of FADC slices summed in extraction - -public: - - MHPedestalPix(const char *name=NULL, const char *title=NULL); - ~MHPedestalPix() {} - - // Setters - void SetNSlices( const Float_t n) { fNSlices = n ; } - - // Getters - Float_t GetNSlices() const { return fNSlices; } - - // Others - void Renorm(); - - ClassDef(MHPedestalPix, 1) // Histogram class for Charge Pedestal Pixel -}; - -#endif Index: trunk/MagicSoft/Mars/manalysis/MHadronness.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHadronness.cc (revision 6973) +++ (revision ) @@ -1,71 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 5/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MHadronness -// -// The Value must be in the range [0,1] -// It should be the value used for gamma/hadron seperation. -// For quality histograms look at MHHadronness -// -///////////////////////////////////////////////////////////////////////////// -#include "MHadronness.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHadronness); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MHadronness::MHadronness(const char *name, const char *title) - : fHadronness(-1) -{ - fName = name ? name : "MHadronness"; - fTitle = title ? title : "Storage container for hadroness value"; -} - -// -------------------------------------------------------------------------- -// -// Reset hadroness, -1 indicates: invalid value -// -void MHadronness::Reset() -{ - fHadronness = -1; -} - -// -------------------------------------------------------------------------- -// -// Print the hillas Parameters to *fLog -// -void MHadronness::Print(Option_t *) const -{ - *fLog << all << GetDescriptor() << ": Hadronness = " << fHadronness << endl; -} - Index: trunk/MagicSoft/Mars/manalysis/MHadronness.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHadronness.h (revision 6973) +++ (revision ) @@ -1,26 +1,0 @@ -#ifndef MARS_MHadronness -#define MARS_MHadronness - -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif - -class MHadronness : public MParContainer -{ -private: - Double_t fHadronness; // [0,1] Hadronness - -public: - MHadronness(const char *name=NULL, const char *title=NULL); - - void Reset(); - - virtual void Print(Option_t *opt=NULL) const; - - Double_t GetHadronness() const { return fHadronness; } - void SetHadronness(Double_t h) { fHadronness = h; } - - ClassDef(MHadronness, 1) // Storage Container for the hadroness -}; - -#endif Index: trunk/MagicSoft/Mars/manalysis/MPedCalUpdate.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedCalUpdate.cc (revision 6973) +++ (revision ) @@ -1,236 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 03/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MPedCalUpdate -// -// This task reads files stored by MJ* storing pedestals (F2) and -// calibration constants (F1) while an eventloop is running. -// -// The update is done according to the run-number. The information -// which ped/cal files corresponds to which run-number is stored in -// an ascii file. In a next step this ascii file can be created -// automatically (hopefully from the Database) -// -///////////////////////////////////////////////////////////////////////////// -#include "MPedCalUpdate.h" - -#include -#include - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" - -#include "MRawRunHeader.h" -#include "MPedPhotCam.h" -#include "MCalibrationChargeCam.h" -#include "MBadPixelsCam.h" - -ClassImp(MPedCalUpdate); - -using namespace std; - -static const TString gsDefName = "MPedCalUpdate"; -static const TString gsDefTitle = ""; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MPedCalUpdate::MPedCalUpdate(const char *fname, const char *name, const char *title) - : fPedPhot(0), fCalCam(0), fFileName(fname) -{ - fName = name ? name : gsDefName.Data(); - fTitle = title ? title : gsDefTitle.Data(); - - fLast=0; -} - -Bool_t MPedCalUpdate::ReadPed(const char *fname) -{ - TFile file(fname, "READ"); - if (fPedPhot->Read()<=0) - { - *fLog << err << "Unable to read MPedPhotCam from " << fname << endl; - return kFALSE; - } - - if (file.FindKey("MBadPixelsCam")) - { - MBadPixelsCam bad; - if (bad.Read()<=0) - { - *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl; - return kFALSE; - } - fBadPix->Merge(bad); - } - return kTRUE; -} - -Bool_t MPedCalUpdate::ReadCal(const char *fname) -{ - TFile file(fname, "READ"); - if (fCalCam->Read()<=0) - { - *fLog << err << "Unable to read MCalibrationChargeCam from " << fname << endl; - return kFALSE; - } - - if (file.FindKey("MBadPixelsCam")) - { - MBadPixelsCam bad; - if (bad.Read()<=0) - { - *fLog << "Unable to read MBadPixelsCam from " << fname << endl; - return kFALSE; - } - fBadPix->Merge(bad); - } - return kTRUE; -} - -Bool_t MPedCalUpdate::ReadPC(TString &line) -{ - const Int_t f = line.First(' '); - if (f<0) - { - *fLog << err << "Entry for pedestal file not found!" << endl; - return kFALSE; - } - - const Int_t l = line.Last(' ')+1; - if (l<0) - { - *fLog << err << "Entry for calibration file not found!" << endl; - return kFALSE; - } - - const TString fname1 = line(0, f); - const TString fname2 = line(l, line.Length()); - - TString mid(line(f, l-f)); - mid=mid.Strip(TString::kBoth); - if (!mid.IsNull()) - { - *fLog << err << "Found three filenames '" << fname1 << "' '" << mid << "' '" << fname2 << "'" << endl; - return kFALSE; - } - - if (!ReadPed(fname1)) - return kFALSE; - - if (!ReadCal(fname2)) - return kFALSE; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// -Bool_t MPedCalUpdate::ReInit(MParList *pList) -{ - MRawRunHeader *header = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!header) - { - *fLog << err << "MRawRunHeader not found... aborting." << endl; - return kFALSE; - } - fPedPhot = (MPedPhotCam*)pList->FindCreateObj("MPedPhotCam"); - if (!fPedPhot) - return kFALSE; - fBadPix = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam"); - if (!fBadPix) - return kFALSE; - fCalCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam"); - if (!fCalCam) - return kFALSE; - - const Int_t run = header->GetRunNumber(); - - ifstream fin(fFileName); - if (!fin) - { - *fLog << err << "Cannot open file: " << fFileName << endl; - return kFALSE; - } - - Int_t last = 0; - - TString readline; - - while (1) - { - // Read line from file - TString line; - line.ReadLine(fin); - if (!fin && !readline.IsNull() && fLast!=last) - { - *fLog << dbg << "Reading line " << readline << endl; - fLast = last; - return ReadPC(readline); - } - - if (!fin) - break; - - line = line.Strip(TString::kBoth); - - // Interprete line - Int_t len, num; - if (sscanf(line.Data(), "%d %n", &num, &len)!=1) - continue; - - // Check whether entries are in increasing order - if (num<=last) - { - *fLog << err << "Runs must be in increasing order..." << endl; - return kFALSE; - } - last = num; - - // Check if a new set of files must be read - if (num<=run) - { - readline=line; - readline.Remove(0, len); - continue; - } - - // This is the case if the correct files were already read - if (fLast==num) - return kTRUE; - - // Read files given in new line - *fLog << dbg << "Reading line " << readline << endl; - fLast = num; - return ReadPC(readline); - } - - return kTRUE; -} Index: trunk/MagicSoft/Mars/manalysis/MPedCalUpdate.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedCalUpdate.h (revision 6973) +++ (revision ) @@ -1,35 +1,0 @@ -#ifndef MARS_MPedCalUpdate -#define MARS_MPedCalUpdate - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class MPedPhotCam; -class MCalibrationChargeCam; -class MBadPixelsCam; - -class MPedCalUpdate : public MTask -{ -private: - MPedPhotCam *fPedPhot; - MCalibrationChargeCam *fCalCam; - MBadPixelsCam *fBadPix; - - TString fFileName; - - Int_t fLast; - - Bool_t ReadCal(const char *name); - Bool_t ReadPed(const char *name); - Bool_t ReadPC(TString &line); - Bool_t ReInit(MParList *pList); - -public: - MPedCalUpdate(const char *filename, const char *name=NULL, const char *title=NULL); - - ClassDef(MPedCalUpdate, 1) // Task to update Pedestals and Cal-Constants while an eventloop is running -}; - -#endif - Index: trunk/MagicSoft/Mars/manalysis/MPedestalCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalCalc.cc (revision 6973) +++ (revision ) @@ -1,229 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Josep Flix 11/2002 -! Author(s): Thomas Bretz 11/2002 -! -! Copyright: MAGIC Software Development, 2002 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MPedestalCalc -// -// Task to calculate the pedestals from the event stream. -// -// From the events two streams of pedestals are created like in the -// following table -// MRawEvtData: 1 2 3 4 5 6 7 8 9 10 11 12 -// MPedestalCam;1: ------1------ ------2------ ------3------... -// MPedestalCam;2: ------1------ ------2------ ------... -// -// Input Containers: -// MRawEvtData -// -// Output Containers: -// MPedestalCam;1 -// MPedestalCam;2 -// -///////////////////////////////////////////////////////////////////////////// -#include "MPedestalCalc.h" - -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MTime.h" -#include "MHFadcCam.h" -#include "MRawEvtPixelIter.h" - -#include "MPedestalCam.h" -#include "MPedestalPix.h" - -ClassImp(MPedestalCalc); - -MPedestalCalc::MPedestalCalc(const char *name, const char *title) - : fNumBranches(2), fTimeSlot(6), fHists(NULL) -{ - fName = name ? name : "MPedestalCalc"; - fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data"; - - AddToBranchList("fHiGainPixId"); - AddToBranchList("fLoGainPixId"); - AddToBranchList("fHiGainFadcSamples"); - AddToBranchList("fLoGainFadcSamples"); -} - -Int_t MPedestalCalc::PreProcess(MParList *pList) -{ - if (fHists) - { - *fLog << err << "ERROR - Previous PostProcess not called." << endl; - return kFALSE; - } - - fHists = new MHFadcCam[fNumBranches]; - fStart = new MTime [fNumBranches]; - - fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); - if (!fRawEvt) - { - *fLog << dbginf << "MRawEvtData not found... aborting." << endl; - return kFALSE; - } - - fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam"); - if (!fPedestals) - return kFALSE; - - fPedTime = (MTime*)pList->FindCreateObj("MTime", "MPedestalTime"); - if (!fPedTime) - return kFALSE; - - // - // Fixme: FindCreateObj --> FindObject - // - fEvtTime = (MTime*)pList->FindCreateObj("MTime", "MRawEvtTime"); - if (!fEvtTime) - return kFALSE; - - for (int i=0; iSetDuration(fTimeSlot); - fEvtTime->SetTime(0, 0); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Checks whether the current event exceed the validity range of a pedestal -// in one of the streams. If this is the case the data from the histogram -// is dumped as pedestals into the corresponding containers and the -// histograms are reset. -// -// Then the current event is filled into the histograms. -// -Int_t MPedestalCalc::Process() -{ - // - // Time when filling of stream a/b must be restarted - // - for (int i=0; iSetTime(fEvtTime->GetTimeLo()+10, 0); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the event in the histogram if the current time of the -// event is _after_ the start time for filling -// -void MPedestalCalc::Fill(Int_t i) -{ - if (*fEvtTime < fStart[i]) - return; - - fHists[i].Fill(fRawEvt); -} - -// -------------------------------------------------------------------------- -// -// If the current time of the event is equal or after the time -// which is given by the start time plus the validity range -// calculate the corresponding set of pedestals and reset the -// histograms -// -void MPedestalCalc::Check(Int_t i) -{ - if (*fEvtTime-fStart[i] < fTimeSlot*10) - return; - - Calc(i); - - fStart[i] = *fEvtTime; - fHists[i].ResetHistograms(); -} - -// -------------------------------------------------------------------------- -// -// Deletes all dynamicly allocated arrays -// -Bool_t MPedestalCalc::PostProcess() -{ - delete fHists; - delete fStart; - - fHists=NULL; - - return kTRUE; -} - -/* - #include "MRawEvtData.h" - */ -// -------------------------------------------------------------------------- -// -// Dumps the mean, rms and their errors from the filled histograms into -// the Pedestal container. Sets the ReadyToSaveFlag of the MPedestalCam -// container. -// -void MPedestalCalc::Calc(Int_t i) const -{ - // FIXME! Fit +- 3 sigma to avoid outliers... - - MRawEvtPixelIter pixel(fRawEvt); - while (pixel.Next()) - { - const UInt_t pixid = pixel.GetPixelId(); - - const TH1 &h = *fHists[i].GetHistHi(pixid); - - const Int_t entries = (Int_t)h.GetEntries(); - const Float_t meanhi = h.GetMean(); - const Float_t rmshi = h.GetRMS(); - - const Float_t meanhierr = rmshi/sqrt(entries); - const Float_t rmshierr = rmshi/sqrt(entries*2); - - MPedestalPix &pix = (*fPedestals)[pixid]; - pix.SetPedestal(meanhi, rmshi); - pix.SetPedestalRms(meanhierr, rmshierr); - } - - *fPedTime = fStart[i]; - - fPedTime->SetReadyToSave(); - fPedestals->SetReadyToSave(); - - /* - *fLog << i << "[" << fHists[i].GetHistHi(0)->GetEntries()/fRawEvt->GetNumHiGainSamples() << "]: "; - *fLog << fEvtTime->GetTimeLo() << ": Calc ["; - *fLog << fStart[i].GetTimeLo() << " < "; - *fLog << fStart[i].GetTimeLo()+fTimeSlot*10 << "]" << endl; - */ -} Index: trunk/MagicSoft/Mars/manalysis/MPedestalCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalCalc.h (revision 6973) +++ (revision ) @@ -1,59 +1,0 @@ -#ifndef MARS_MPedestalCalc -#define MARS_MPedestalCalc - -///////////////////////////////////////////////////////////////////////////// -// // -// MPedestalCalc // -// // -// Implementation of ped. Eval. defined in Jan 02 // -// // -///////////////////////////////////////////////////////////////////////////// -#ifndef MARS_MTask -#include "MTask.h" -#endif - -#ifndef MARS_MHFadCam -#include "MHFadcCam.h" -#endif - -class MTime; -class MRawEvtData; -class MPedestalCam; - -class MPedestalCalc : public MTask -{ -private: - Int_t fNumBranches; - Float_t fTimeSlot; - - MRawEvtData *fRawEvt; - MPedestalCam *fPedestals; // - - MHFadcCam *fHists; //[fNumBranches] - - MTime *fPedTime; // time of the pedestal event - MTime *fEvtTime; // time of the current event - MTime *fStart; //[fNumBranches] starting time of the current pedestal - - - void Calc(Int_t i) const; - void Fill(Int_t i); - void Check(Int_t i); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - -public: - MPedestalCalc(const char *name=NULL, const char *title=NULL); - - void SetTimeSlot(Float_t newslot) { fTimeSlot = newslot; } - void SetNumBranches(Int_t num) { fNumBranches = num; } - - ClassDef(MPedestalCalc, 0) // Task to calculate the pestels from pedestal events -}; - -#endif - - - Index: trunk/MagicSoft/Mars/manalysis/MPedestalSum.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalSum.cc (revision 6973) +++ (revision ) @@ -1,227 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Abelardo Moralejo 7/2002 -! Author(s): Thomas Bretz 2002 -! -! Copyright: MAGIC Software Development, 2002-2003 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MPedestalSum -// -// This is a task which calculates the number of photons from the FADC -// time slices. It weights the each slice according to the numbers in -// the array fWeight (default: all slices added up with weight 1). -// -// The weights are rescaled, such that sum(weigths)=num slices -// -// Input Containers: -// MRawRunHeader, MRawEvtData, MPedestalCam -// -// Output Containers: -// MCerPhotEvt -// -////////////////////////////////////////////////////////////////////////////// -#include "MPedestalSum.h" - -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MMcRunHeader.hxx" - -#include "MRawRunHeader.h" -#include "MRawEvtData.h" // MRawEvtData::GetNumPixels -#include "MCerPhotEvt.h" -#include "MPedestalPix.h" -#include "MPedestalCam.h" -#include "MRawEvtPixelIter.h" - -ClassImp(MPedestalSum); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MPedestalSum::MPedestalSum(const char *name, const char *title) -{ - fName = name ? name : "MPedestalSum"; - fTitle = title ? title : "Calculate pixel signal from FADC data"; - - AddToBranchList("MRawEvtData.fHiGainPixId"); - AddToBranchList("MRawEvtData.fLoGainPixId"); - AddToBranchList("MRawEvtData.fHiGainFadcSamples"); - AddToBranchList("MRawEvtData.fLoGainFadcSamples"); - - SetDefaultWeights(); -} - -// -------------------------------------------------------------------------- -// -// The PreProcess searches for the following input containers: -// - MRawRunHeader -// - MRawEvtData -// - MPedestalCam -// -// The following output containers are also searched and created if -// they were not found: -// - MCerPhotEvt -// -Int_t MPedestalSum::PreProcess(MParList *pList) -{ - fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!fRunHeader) - { - *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; - return kFALSE; - } - fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); - if (!fRawEvt) - { - *fLog << dbginf << "MRawEvtData not found... aborting." << endl; - return kFALSE; - } - - fCerPhotEvt = (MCerPhotEvt*)pList->FindCreateObj("MCerPhotEvt"); - if (!fCerPhotEvt) - { - *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; - return kFALSE; - } - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Check for the run type and camera version. -// If the file is a MC file and the used camera version is <= 40 -// we enable a fix for truncated pedestal means in this version. -// -Bool_t MPedestalSum::ReInit(MParList *pList) -{ - /* - const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!runheader) - { - *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; - return kTRUE; - } - - if (fRunHeader->GetNumSamplesHiGain() != fWeight.GetSize()) - { - *fLog << dbginf << "Number of FADC slices (" << fRunHeader->GetNumSamplesHiGain() <<") is different from assumed one (" << fWeight.GetSize() << ")... aborting." << endl; - return kFALSE; - } - - Bool_t fIsMcFile = runheader->GetRunType() == kRTMonteCarlo; - if (!fIsMcFile) - return kTRUE; - - ScalePedestals(); - - MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader"); - if (!mcrunheader) - { - *fLog << warn << dbginf << "Warning - cannot check for camera version, MC run header not found." << endl; - return kTRUE; - } - - fEnableFix = kFALSE; - if (mcrunheader->GetCamVersion() <= 40) - fEnableFix = kTRUE; - */ - return kTRUE; -} - -void MPedestalSum::ScalePedestals() -{ - Int_t n = 577; //FIXME: fPedestals->GetNumPixel(); - - for (int idx=0; idxSetReadyToSave(); -} - -// -------------------------------------------------------------------------- -// -// Calculate the integral of the FADC time slices and store them as a new -// pixel in the MCerPhotEvt container. -// -Int_t MPedestalSum::Process() -{ - //fCerPhotEvt->InitSize(fRawEvt->GetNumPixels()); - -// if (fIsMcFile) -// ScalePedestals(); - - MRawEvtPixelIter pixel(fRawEvt); - - while (pixel.Next()) - { - const UInt_t idx = pixel.GetPixelId(); - - const int n = fRunHeader->GetNumSamplesHiGain(); - - // - // Calculate pixel signal unless it has all FADC slices empty: - // - Byte_t *ptr = pixel.GetHiGainSamples(); - - Float_t nphot = 0; - Float_t nerr = 0; - for(Int_t i=0; iAddPixel(idx, nphot, sqrt(nerr/n-nphot*nphot)); - - // FIXME! Handling of Lo Gains is missing! - } - - fCerPhotEvt->FixSize(); - fCerPhotEvt->SetReadyToSave(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Set default values for the number of slices and weights: -// -void MPedestalSum::SetDefaultWeights() -{ - const Float_t dummy[15] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }; - fWeight.Set(15, dummy); -} Index: trunk/MagicSoft/Mars/manalysis/MPedestalSum.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalSum.h (revision 6973) +++ (revision ) @@ -1,57 +1,0 @@ -#ifndef MARS_MPedestalSum -#define MARS_MPedestalSum - -///////////////////////////////////////////////////////////////////////////// -// // -// MPedestalSum // -// // -// Integrates the desired ADC time slices of one pixel and substracts the // -// pedestal (offset) value // -// // -///////////////////////////////////////////////////////////////////////////// -#ifndef ROOT_TArrayF -#include -#endif - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class MRawEvtData; -class MPedestalCam; -class MCerPhotEvt; -class MRawRunHeader; - -class MPedestalSum : public MTask -{ - MPedestalCam *fPedestals; // Pedestals of all pixels in the camera - MRawEvtData *fRawEvt; // raw event data (time slices) - MCerPhotEvt *fCerPhotEvt; // Cerenkov Photon Event used for calculation - MRawRunHeader *fRunHeader; // RunHeader information - - Bool_t fEnableFix; // fix for a bug in files from older camera versions (<=40) - Bool_t fIsMcFile; - - TArrayF fWeight; // Weights for adding up the ADC slices - Float_t fSumWeights; - Float_t fSumQuadWeights; - - void SetDefaultWeights(); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - - Bool_t ReInit(MParList *pList); - void ScalePedestals(); - -public: - MPedestalSum(const char *name=NULL, const char *title=NULL); - - // FIXME: The array size should be checked! - void SetWeights(const TArrayF &w) { fWeight = w; } - - ClassDef(MPedestalSum, 0) // Task to calculate cerenkov photons from raw data -}; - - -#endif Index: trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc (revision 6973) +++ (revision ) @@ -1,552 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -! Wolfgang Wittek 12/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MPedestalWorkaround // -// // -// copies the pedestal value and the pedestal RMS // -// from the container MPedPhotCam // -// into the contaier MPedestalCam // -// // -// put the zenith angle into MMcEvt; // -// the zenith angle is taken from the runbooks // -// (1 fixed zenith angle for a given run) // -// // -// this workaround is necessary // -// - as long as the analysis classes // -// take the pedestals from MPedestalCam; in the long run they have to // -// be taken from MPedPhotCam // -// - and as long as the container for the zenith angle is not defined // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MPedestalWorkaround.h" - -#include -#include "MLog.h" -#include "MLogManip.h" -#include "MParList.h" -#include "MGeomCam.h" - -#include "MPedPhotCam.h" -#include "MPedPhotPix.h" -#include "MMcEvt.hxx" -#include "MRawRunHeader.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; - } - - fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!fRun) - { - *fLog << err << "MRawRunHeader not found... aborting." << endl; - return kFALSE; - } - - fMcEvt = (MMcEvt*)pList->FindCreateObj("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << "MMcEvt not found... aborting." << endl; - return kFALSE; - } - - return kTRUE; -} - -// ------------------------------------------------------------------------ -// -Int_t MPedestalWorkaround::Process() -{ - //------------------------------------------------------------------- - // copy the pedestal and the pedestal RMS - // from MPedPhotCam into MPedestalCam - - // set pedestalRMS = 0 if it is too high - UInt_t imaxnumpix = fCam->GetNumPixels(); - for (UInt_t i=0; i 10.0) - (*fPedPhot)[i].SetRms(0.0); - } - - - //------------------------------------------------------------------- - // put the zenith angle into MMcEvt - - Double_t thetadeg; - Double_t thetarad; - - - /* - Int_t run = fRun->GetRunNumber(); - - if (run == 3127) thetadeg = 27.2; // Crab - else if (run == 3128) thetadeg = 25.6; - else if (run == 3129) thetadeg = 24.3; - else if (run == 3130) thetadeg = 23.9; - - - else if (run == 3132) thetadeg = 22.9; - else if (run == 3133) thetadeg = 22.8; - else if (run == 3134) thetadeg = 22.3; - else if (run == 3135) thetadeg = 21.9; - else if (run == 3136) thetadeg = 21.5; - else if (run == 3137) thetadeg = 21.1; - else if (run == 3138) thetadeg = 20.8; - else if (run == 3139) thetadeg = 20.4; - - else if (run == 3140) thetadeg = 19.5; - else if (run == 3141) thetadeg = 19.4; - else if (run == 3142) thetadeg = 19.0; - else if (run == 3143) thetadeg = 18.6; - else if (run == 3144) thetadeg = 13.0; - else if (run == 3145) thetadeg = 12.4; - else if (run == 3146) thetadeg = 12.1; - else if (run == 3147) thetadeg = 11.7; - else if (run == 3148) thetadeg = 11.3; - else if (run == 3149) thetadeg = 11.9; - - else if (run == 3150) thetadeg = 10.6; - else if (run == 3151) thetadeg = 10.3; - else if (run == 3152) thetadeg = 10.0; - else if (run == 3153) thetadeg = 9.6; - else if (run == 3154) thetadeg = 9.3; - else if (run == 3155) thetadeg = 9.0; - else if (run == 3156) thetadeg = 8.7; - else if (run == 3157) thetadeg = 8.4; - else if (run == 3158) thetadeg = 8.1; - else if (run == 3159) thetadeg = 7.9; - - else if (run == 3160) thetadeg = 7.7; - else if (run == 3161) thetadeg = 7.3; - else if (run == 3162) thetadeg = 7.2; - else if (run == 3163) thetadeg = 7.0; - else if (run == 3164) thetadeg = 6.8; - else if (run == 3165) thetadeg = 6.7; - else if (run == 3166) thetadeg = 6.6; - else if (run == 3167) thetadeg = 6.5; - else if (run == 3168) thetadeg = 6.4; - else if (run == 3169) thetadeg = 6.4; - - else if (run == 3170) thetadeg = 6.4; - else if (run == 3171) thetadeg = 6.4; - else if (run == 3172) thetadeg = 6.5; - else if (run == 3173) thetadeg = 6.6; - else if (run == 3174) thetadeg = 6.7; - - else if (run == 3176) thetadeg = 7.1; - else if (run == 3177) thetadeg = 7.4; - else if (run == 3178) thetadeg = 7.6; - else if (run == 3179) thetadeg = 7.9; - - - else if (run == 3182) thetadeg = 8.4; - else if (run == 3183) thetadeg = 8.9; - else if (run == 3184) thetadeg = 9.2; - else if (run == 3185) thetadeg = 9.5; - else if (run == 3186) thetadeg = 9.8; - else if (run == 3187) thetadeg = 10.5; - else if (run == 3188) thetadeg = 10.9; - else if (run == 3189) thetadeg = 11.2; - - else if (run == 3190) thetadeg = 11.6; - else if (run == 3191) thetadeg = 11.6; - else if (run == 3192) thetadeg = 12.4; - else if (run == 3193) thetadeg = 12.7; - else if (run == 3194) thetadeg = 13.1; - else if (run == 3195) thetadeg = 13.5; - else if (run == 3196) thetadeg = 13.9; - else if (run == 3197) thetadeg = 14.3; - else if (run == 3198) thetadeg = 14.7; - else if (run == 3199) thetadeg = 15.1; - - else if (run == 3200) thetadeg = 15.6; - else if (run == 3201) thetadeg = 16.0; - else if (run == 3202) thetadeg = 16.5; - else if (run == 3203) thetadeg = 16.9; - else if (run == 3204) thetadeg = 17.3; - else if (run == 3205) thetadeg = 17.7; - else if (run == 3206) thetadeg = 18.2; - else if (run == 3207) thetadeg = 18.6; - else if (run == 3208) thetadeg = 19.0; - else if (run == 3209) thetadeg = 19.4; - - else if (run == 3210) thetadeg = 19.9; - else if (run == 3211) thetadeg = 20.4; - else if (run == 3212) thetadeg = 20.8; - else if (run == 3213) thetadeg = 21.2; - else if (run == 3214) thetadeg = 21.7; - else if (run == 3215) thetadeg = 22.2; - else if (run == 3216) thetadeg = 25.6; // Off Crab1 - else if (run == 3217) thetadeg = 25.0; - else if (run == 3218) thetadeg = 24.5; - else if (run == 3219) thetadeg = 24.0; - - else if (run == 3220) thetadeg = 23.5; - else if (run == 3221) thetadeg = 22.5; - else if (run == 3222) thetadeg = 22.1; - - else if (run == 3225) thetadeg = 15.1; - else if (run == 3226) thetadeg = 15.0; - else if (run == 3227) thetadeg = 14.5; - else if (run == 3228) thetadeg = 14.1; - else if (run == 3229) thetadeg = 13.8; - - else if (run == 3230) thetadeg = 13.3; - else if (run == 3231) thetadeg = 13.0; - else if (run == 3232) thetadeg = 12.6; - else if (run == 3233) thetadeg = 12.3; - else if (run == 3234) thetadeg = 12.0; - else if (run == 3235) thetadeg = 11.6; - else if (run == 3236) thetadeg = 11.3; - else if (run == 3237) thetadeg = 11.0; - else if (run == 3238) thetadeg = 10.8; - else if (run == 3239) thetadeg = 10.4; - - else if (run == 3240) thetadeg = 10.1; - else if (run == 3241) thetadeg = 9.9; - else if (run == 3242) thetadeg = 9.6; - else if (run == 3243) thetadeg = 9.4; - else if (run == 3244) thetadeg = 9.2; - else if (run == 3245) thetadeg = 9.0; - else if (run == 3246) thetadeg = 8.9; - else if (run == 3247) thetadeg = 8.8; - else if (run == 3248) thetadeg = 8.7; - else if (run == 3249) thetadeg = 8.6; - - else if (run == 3250) thetadeg = 8.6; - else if (run == 3251) thetadeg = 8.6; - else if (run == 3252) thetadeg = 8.6; - else if (run == 3253) thetadeg = 8.7; - else if (run == 3254) thetadeg = 8.8; - else if (run == 3255) thetadeg = 8.9; - else if (run == 3256) thetadeg = 9.1; - else if (run == 3257) thetadeg = 9.3; - else if (run == 3258) thetadeg = 9.5; - else if (run == 3259) thetadeg = 9.7; - else if (run == 3260) thetadeg = 9.9; - - else if (run == 3261) thetadeg = 10.2; - else if (run == 3262) thetadeg = 10.5; - else if (run == 3263) thetadeg = 10.8; - else if (run == 3264) thetadeg = 11.1; - else if (run == 3265) thetadeg = 11.4; - else if (run == 3266) thetadeg = 11.8; - else if (run == 3267) thetadeg = 12.1; - else if (run == 3268) thetadeg = 12.5; - else if (run == 3269) thetadeg = 12.8; - - else if (run == 3270) thetadeg = 13.2; - else if (run == 3271) thetadeg = 13.5; - else if (run == 3272) thetadeg = 13.9; - else if (run == 3273) thetadeg = 14.0; - else if (run == 3274) thetadeg = 14.4; - - else if (run == 3284) thetadeg = 7.0; // Crab - else if (run == 3285) thetadeg = 7.2; - else if (run == 3286) thetadeg = 7.4; - else if (run == 3287) thetadeg = 7.5; - else if (run == 3288) thetadeg = 8.4; - else if (run == 3289) thetadeg = 9.0; - - else if (run == 3290) thetadeg = 9.4; - else if (run == 3291) thetadeg = 9.8; - else if (run == 3292) thetadeg = 10.2; - else if (run == 3293) thetadeg = 10.5; - else if (run == 3294) thetadeg = 10.9; - else if (run == 3295) thetadeg = 11.3; - else if (run == 3296) thetadeg = 11.8; - else if (run == 3297) thetadeg = 12.2; - else if (run == 3298) thetadeg = 12.6; - else if (run == 3299) thetadeg = 13.0; - - else if (run == 3300) thetadeg = 13.5; - else if (run == 3301) thetadeg = 13.9; - else if (run == 3302) thetadeg = 14.3; - else if (run == 3303) thetadeg = 14.8; - else if (run == 3304) thetadeg = 15.2; - else if (run == 3305) thetadeg = 15.7; - else if (run == 3306) thetadeg = 16.2; - else if (run == 3307) thetadeg = 16.6; - else if (run == 3308) thetadeg = 17.1; - else if (run == 3309) thetadeg = 17.6; - - else if (run == 3310) thetadeg = 17.9; - else if (run == 3311) thetadeg = 18.4; - else if (run == 3312) thetadeg = 18.9; - else if (run == 3313) thetadeg = 19.3; - else if (run == 3314) thetadeg = 19.8; - else if (run == 3315) thetadeg = 20.1; - else if (run == 3316) thetadeg = 20.7; - else if (run == 3317) thetadeg = 21.2; - else if (run == 3318) thetadeg = 21.7; - else if (run == 3319) thetadeg = 22.1; - - else if (run == 3320) thetadeg = 22.5; - else if (run == 3321) thetadeg = 23.1; - else if (run == 3322) thetadeg = 23.6; - else if (run == 3323) thetadeg = 24.1; - else if (run == 3324) thetadeg = 24.6; - else if (run == 3325) thetadeg = 24.9; - else if (run == 3326) thetadeg = 25.5; - else if (run == 3327) thetadeg = 26.0; - else if (run == 3328) thetadeg = 26.0; - else if (run == 3329) thetadeg = 26.6; - - else if (run == 3330) thetadeg = 26.6; - else if (run == 3331) thetadeg = 27.1; - else if (run == 3332) thetadeg = 27.7; - else if (run == 3333) thetadeg = 28.2; - else if (run == 3334) thetadeg = 28.5; - - else if (run == 3340) thetadeg = 10.5; - else if (run == 3341) thetadeg = 10.3; - else if (run == 3342) thetadeg = 9.6; - else if (run == 3343) thetadeg = 9.2; - else if (run == 3344) thetadeg = 8.9; - else if (run == 3345) thetadeg = 8.6; - else if (run == 3346) thetadeg = 8.3; - else if (run == 3347) thetadeg = 8.0; - else if (run == 3348) thetadeg = 7.7; - else if (run == 3349) thetadeg = 7.5; - - else if (run == 3350) thetadeg = 7.2; - else if (run == 3351) thetadeg = 7.0; - else if (run == 3352) thetadeg = 6.8; - else if (run == 3353) thetadeg = 6.7; - else if (run == 3354) thetadeg = 6.6; - else if (run == 3355) thetadeg = 6.5; - else if (run == 3356) thetadeg = 6.4; - else if (run == 3357) thetadeg = 6.4; - else if (run == 3358) thetadeg = 6.4; - else if (run == 3359) thetadeg = 6.5; - - else if (run == 3360) thetadeg = 6.6; - - else if (run == 3362) thetadeg = 6.7; - else if (run == 3363) thetadeg = 6.8; - else if (run == 3364) thetadeg = 7.0; - else if (run == 3365) thetadeg = 7.2; - else if (run == 3366) thetadeg = 7.5; - else if (run == 3367) thetadeg = 7.7; - else if (run == 3368) thetadeg = 8.0; - else if (run == 3369) thetadeg = 8.3; - - else if (run == 3370) thetadeg = 8.6; - else if (run == 3371) thetadeg = 9.0; - else if (run == 3372) thetadeg = 9.3; - else if (run == 3373) thetadeg = 9.6; - else if (run == 3374) thetadeg = 10.0; - else if (run == 3375) thetadeg = 10.4; - else if (run == 3376) thetadeg = 10.7; - else if (run == 3377) thetadeg = 11.1; - else if (run == 3378) thetadeg = 11.5; - else if (run == 3379) thetadeg = 11.9; - - else if (run == 3380) thetadeg = 12.3; - else if (run == 3381) thetadeg = 12.7; - else if (run == 3382) thetadeg = 13.1; - else if (run == 3383) thetadeg = 13.5; - else if (run == 3384) thetadeg = 13.9; - else if (run == 3385) thetadeg = 14.3; - else if (run == 3386) thetadeg = 14.7; - else if (run == 3387) thetadeg = 15.2; - else if (run == 3388) thetadeg = 15.6; - else if (run == 3389) thetadeg = 16.0; - - else if (run == 3390) thetadeg = 16.4; - else if (run == 3391) thetadeg = 16.7; - else if (run == 3392) thetadeg = 17.9; - else if (run == 3393) thetadeg = 18.3; - else if (run == 3394) thetadeg = 18.7; - else if (run == 3395) thetadeg = 19.2; - else if (run == 3396) thetadeg = 19.6; - else if (run == 3397) thetadeg = 20.0; - else if (run == 3398) thetadeg = 20.5; - else if (run == 3399) thetadeg = 20.9; - - else if (run == 3400) thetadeg = 21.4; - else if (run == 3401) thetadeg = 21.8; - else if (run == 3402) thetadeg = 22.1; - else if (run == 3403) thetadeg = 22.6; - else if (run == 3404) thetadeg = 23.1; - else if (run == 3405) thetadeg = 23.4; - else if (run == 3406) thetadeg = 23.9; - else if (run == 3407) thetadeg = 24.3; - else if (run == 3408) thetadeg = 24.7; - else if (run == 3409) thetadeg = 26.9; - - else if (run == 3410) thetadeg = 27.3; - else if (run == 3411) thetadeg = 27.7; - else if (run == 3412) thetadeg = 28.2; - else if (run == 3413) thetadeg = 28.7; - else if (run == 3414) thetadeg = 29.1; - else if (run == 3415) thetadeg = 29.2; - else if (run == 3416) thetadeg = 30.0; - else if (run == 3417) thetadeg = 18.5; // Off Crab1 - else if (run == 3418) thetadeg = 18.4; - else if (run == 3419) thetadeg = 17.5; - - else if (run == 3420) thetadeg = 17.2; - else if (run == 3421) thetadeg = 16.8; - else if (run == 3422) thetadeg = 16.4; - else if (run == 3423) thetadeg = 16.0; - else if (run == 3424) thetadeg = 15.6; - else if (run == 3425) thetadeg = 15.3; - else if (run == 3426) thetadeg = 14.9; - else if (run == 3427) thetadeg = 14.5; - else if (run == 3428) thetadeg = 14.1; - else if (run == 3429) thetadeg = 13.7; - - else if (run == 3430) thetadeg = 13.4; - else if (run == 3431) thetadeg = 13.0; - else if (run == 3432) thetadeg = 12.7; - else if (run == 3433) thetadeg = 12.3; - else if (run == 3434) thetadeg = 12.0; - else if (run == 3435) thetadeg = 12.0; - else if (run == 3436) thetadeg = 11.6; - else if (run == 3437) thetadeg = 11.3; - else if (run == 3438) thetadeg = 11.0; - else if (run == 3439) thetadeg = 10.7; - - else if (run == 3440) thetadeg = 10.4; - else if (run == 3441) thetadeg = 10.1; - else if (run == 3442) thetadeg = 9.9; - else if (run == 3443) thetadeg = 9.8; - else if (run == 3444) thetadeg = 30.8; // Mkn 421 - else if (run == 3445) thetadeg = 30.6; - else if (run == 3446) thetadeg = 29.7; - else if (run == 3447) thetadeg = 29.3; - else if (run == 3448) thetadeg = 28.9; - else if (run == 3449) thetadeg = 28.5; - - else if (run == 3450) thetadeg = 28.1; - else if (run == 3451) thetadeg = 27.7; - else if (run == 3452) thetadeg = 27.3; - else if (run == 3453) thetadeg = 26.9; - else if (run == 3454) thetadeg = 26.5; - else if (run == 3455) thetadeg = 26.1; - else if (run == 3456) thetadeg = 25.7; - else if (run == 3457) thetadeg = 25.3; - else if (run == 3458) thetadeg = 24.9; - else if (run == 3459) thetadeg = 24.5; - - else if (run == 3460) thetadeg = 24.1; - - else if (run == 3463) thetadeg = 23.2; - else if (run == 3464) thetadeg = 22.8; - else if (run == 3465) thetadeg = 22.4; - else if (run == 3466) thetadeg = 22.0; - else if (run == 3467) thetadeg = 21.6; - else if (run == 3468) thetadeg = 21.2; - else if (run == 3469) thetadeg = 20.8; - - else if (run == 3470) thetadeg = 20.4; - else if (run == 3471) thetadeg = 20.1; - else if (run == 3472) thetadeg = 19.7; - else if (run == 3473) thetadeg = 19.3; - else if (run == 3474) thetadeg = 18.9; - else if (run == 3475) thetadeg = 18.5; - else if (run == 3476) thetadeg = 18.2; - else if (run == 3477) thetadeg = 18.1; - else if (run == 3478) thetadeg = 17.7; - - - else if (run == 3480) thetadeg = 17.5; - else if (run == 3481) thetadeg = 17.1; - else if (run == 3482) thetadeg = 16.7; - else if (run == 3483) thetadeg = 16.3; - else if (run == 3484) thetadeg = 16.0; - else if (run == 3485) thetadeg = 15.6; - else if (run == 3486) thetadeg = 15.3; - else if (run == 3487) thetadeg = 15.0; - else if (run == 3488) thetadeg = 14.8; - else if (run == 3489) thetadeg = 14.8; - - else - { - *fLog << warn << "MPedestalWorkaround : no zenith angle for run number = " - << run << endl; - thetadeg = 90.0; - } - - thetarad = thetadeg / kRad2Deg; - */ - - thetadeg = 10.0; - thetarad = thetadeg / kRad2Deg; - fMcEvt->SetTelescopeTheta(thetarad); - - return kTRUE; -} - -Int_t MPedestalWorkaround::PostProcess() -{ - return kTRUE; -} - - - - - - - Index: trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h (revision 6973) +++ (revision ) @@ -1,47 +1,0 @@ -#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 MMcEvt; -class MRawRunHeader; - -class MPedestalWorkaround : public MTask -{ -private: - - MPedestalCam *fPed; - MPedPhotCam *fPedPhot; - MGeomCam *fCam; - MMcEvt *fMcEvt; - MRawRunHeader *fRun; - - 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/MPointingCorr.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPointingCorr.cc (revision 6973) +++ (revision ) @@ -1,134 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 03/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MPointingCorr // -// // -// This is a task to do the pointing correction // -// // -///////////////////////////////////////////////////////////////////////////// - -#include "MPointingCorr.h" - -#include "MParList.h" - -#include "MMcEvt.hxx" -#include "MSrcPosCam.h" -#include "MGeomCam.h" -#include "MParameters.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MPointingCorr); - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MPointingCorr::MPointingCorr(const char *name, const char *title, - const char *srcname) -{ - fName = name ? name : "MPointingCorr"; - fTitle = title ? title : "Task to do the pointing correction"; - - fSrcName = srcname; -} - -// -------------------------------------------------------------------------- -// -// -// -// -// -Int_t MPointingCorr::PreProcess(MParList *pList) -{ - MGeomCam *geom = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!geom) - *fLog << warn << GetDescriptor() << ": No Camera Geometry available. Using mm-scale for histograms." << endl; - else - { - fMm2Deg = geom->GetConvMm2Deg(); - } - - - fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << dbginf << "MMcEvt not found... aborting." << endl; - return kFALSE; - } - - fHourAngle = (MParameterD*)pList->FindObject("HourAngle", "MParameterD"); - if (!fHourAngle) - { - *fLog << dbginf << "HourAngle not found... aborting." << endl; - return kFALSE; - } - - - fSrcPos = (MSrcPosCam*)pList->FindObject(fSrcName, "MSrcPosCam"); - if (!fSrcPos) - { - *fLog << err << dbginf << fSrcName << " [MSrcPosCam] not found... aborting." << endl; - return kFALSE; - } - - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Do the pointing correction -// -// the parametrization is for Mkn421 2001 data (Daniel Kranich) -// -Int_t MPointingCorr::Process() -{ - // fhourangle is the hour angle [degrees] - // (cx, cy) is the source position in the camera [mm] - // - Float_t fhourangle = fHourAngle->GetVal(); - - //*fLog << "MPointingCorr::Process; fhourangle = " - // << fhourangle << endl; - - Float_t cx = -0.05132 - 0.001064 * fhourangle - - 3.530e-6 * fhourangle * fhourangle; - cx /= fMm2Deg; - - Float_t cy = -0.04883 - 0.0003175* fhourangle - - 2.165e-5 * fhourangle * fhourangle; - cy /= fMm2Deg; - - fSrcPos->SetXY(cx, cy); - - //*fLog << "MPointingCorr::Process; fhourangle, cx, cy, fMm2Deg = " - // << fhourangle << ", " << cx << ", " << cy << ", " - // << fMm2Deg << endl; - - return kTRUE; -} Index: trunk/MagicSoft/Mars/manalysis/MPointingCorr.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPointingCorr.h (revision 6973) +++ (revision ) @@ -1,40 +1,0 @@ -#ifndef MARS_MPointingCorr -#define MARS_MPointingCorr - -///////////////////////////////////////////////////////////////////////////// -// // -// MPointingCorr // -// // -// Task to do the pointing correction // -// // -///////////////////////////////////////////////////////////////////////////// - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class MMcEvt; -class MSrcPosCam; -class MParameterD; - -class MPointingCorr : public MTask -{ -private: - MMcEvt *fMcEvt; - MSrcPosCam *fSrcPos; - TString fSrcName; - MParameterD *fHourAngle; - - Float_t fMm2Deg; - - Int_t PreProcess(MParList *pList); - Int_t Process(); - -public: - MPointingCorr(const char *name=NULL, const char *title=NULL, - const char *srcname="MSrcPosCam"); - - ClassDef(MPointingCorr, 0) // Task to do the pointing correction -}; - -#endif Index: trunk/MagicSoft/Mars/manalysis/MSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSupercuts.cc (revision 6973) +++ (revision ) @@ -1,388 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 the default values - - fLengthUp[0] = 0.2; - fLengthUp[1] = 0.0; - fLengthUp[2] = 0.0; - fLengthUp[3] = 0.0; - fLengthUp[4] = 0.0; - fLengthUp[5] = 0.0; - fLengthUp[6] = 0.0; - fLengthUp[7] = 0.0; - - fLengthLo[0] = 0.; - fLengthLo[1] = 0.; - fLengthLo[2] = 0.; - fLengthLo[3] = 0.; - fLengthLo[4] = 0.; - fLengthLo[5] = 0.; - fLengthLo[6] = 0.; - fLengthLo[7] = 0.; - - fWidthUp[0] = 0.1; - fWidthUp[1] = 0.0; - fWidthUp[2] = 0.0; - fWidthUp[3] = 0.0; - fWidthUp[4] = 0.0; - fWidthUp[5] = 0.0; - fWidthUp[6] = 0.0; - fWidthUp[7] = 0.0; - - fWidthLo[0] = 0.; - fWidthLo[1] = 0.; - fWidthLo[2] = 0.; - fWidthLo[3] = 0.; - fWidthLo[4] = 0.; - fWidthLo[5] = 0.; - fWidthLo[6] = 0.; - fWidthLo[7] = 0.; - - fDistUp[0] = 1.e10; - fDistUp[1] = 0.0; - fDistUp[2] = 0.0; - fDistUp[3] = 0.0; - fDistUp[4] = 0.0; - fDistUp[5] = 0.0; - fDistUp[6] = 0.0; - fDistUp[7] = 0.0; - - fDistLo[0] = 0.0; - fDistLo[1] = 0.0; - fDistLo[2] = 0.0; - fDistLo[3] = 0.0; - fDistLo[4] = 0.0; - fDistLo[5] = 0.0; - fDistLo[6] = 0.0; - fDistLo[7] = 0.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.02; - fStepsizes[1] = 0.0; - fStepsizes[2] = 0.0; - fStepsizes[3] = 0.0; - fStepsizes[4] = 0.0; - fStepsizes[5] = 0.0; - fStepsizes[6] = 0.0; - fStepsizes[7] = 0.0; - - // LengthLo - fStepsizes[8] = 0.0; - fStepsizes[9] = 0.0; - fStepsizes[10] = 0.0; - fStepsizes[11] = 0.0; - fStepsizes[12] = 0.0; - fStepsizes[13] = 0.0; - fStepsizes[14] = 0.0; - fStepsizes[15] = 0.0; - - // WidthUp - fStepsizes[16] = 0.01; - fStepsizes[17] = 0.0; - fStepsizes[18] = 0.0; - fStepsizes[19] = 0.0; - fStepsizes[20] = 0.0; - fStepsizes[21] = 0.0; - fStepsizes[22] = 0.0; - fStepsizes[23] = 0.0; - - // WidthLo - fStepsizes[24] = 0.0; - fStepsizes[25] = 0.0; - fStepsizes[26] = 0.0; - fStepsizes[27] = 0.0; - fStepsizes[28] = 0.0; - fStepsizes[29] = 0.0; - fStepsizes[30] = 0.0; - fStepsizes[31] = 0.0; - - // DistUp - fStepsizes[32] = 0.0; - fStepsizes[33] = 0.0; - fStepsizes[34] = 0.0; - fStepsizes[35] = 0.0; - fStepsizes[36] = 0.0; - fStepsizes[37] = 0.0; - fStepsizes[38] = 0.0; - fStepsizes[39] = 0.0; - - // DistLo - fStepsizes[40] = 0.0; - fStepsizes[41] = 0.0; - fStepsizes[42] = 0.0; - fStepsizes[43] = 0.0; - fStepsizes[44] = 0.0; - fStepsizes[45] = 0.0; - fStepsizes[46] = 0.0; - 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 6973) +++ (revision ) @@ -1,69 +1,0 @@ -#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/mfilter/MFAlpha.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFAlpha.cc (revision 6973) +++ (revision ) @@ -1,129 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 01/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MFAlpha // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MFAlpha.h" - -#include -#include - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" - -#include "MHillasSrc.h" - -ClassImp(MFAlpha); - -using namespace std; - -// -------------------------------------------------------------------------- -// -MFAlpha::MFAlpha(const char *cname, const char type, const Float_t val, - const char *name, const char *title) : fHillas(NULL) -{ - fContName = cname; - Init(type, val, name, title); -} - -// -------------------------------------------------------------------------- -// -MFAlpha::MFAlpha(MHillasSrc *hillas, const char type, const Float_t val, - const char *name, const char *title) : fHillas(hillas) -{ - Init(type, val, name, title); -} - -// -------------------------------------------------------------------------- -// -void MFAlpha::Init(const char type, const Float_t val, - const char *name, const char *title) -{ - fName = name ? name : "MFAlpha"; - fTitle = title ? title : "Filter using the alpha angle"; - - fFilterType = (type=='<' ? kELowerThan : kEGreaterThan); - - if (type!='<' && type!='>') - *fLog << warn << dbginf << "Warning: Neither '<' nor '>' specified... using '>'." << endl; - - fValue = val; - - AddToBranchList(Form("%s.fAlpha", (const char*)fContName)); -} - -// -------------------------------------------------------------------------- -// -Int_t MFAlpha::PreProcess(MParList *pList) -{ - if (fHillas) - return kTRUE; - - fHillas = (MHillasSrc*)pList->FindObject(fContName); - if (fHillas) - return kTRUE; - - *fLog << err << dbginf << fContName << " [MHillas] not found... aborting." << endl; - return kFALSE; -} - -// -------------------------------------------------------------------------- -// -Int_t MFAlpha::Process() -{ - const Float_t alpha = fabs(fHillas->GetAlpha()); - - switch (fFilterType) - { - case kELowerThan: - fResult = (alpha < fValue); - break; - case kEGreaterThan: - fResult = (alpha > fValue); - break; - } - - return kTRUE; -} - -void MFAlpha::StreamPrimitive(ofstream &out) const -{ - if (fHillas) - fHillas->SavePrimitive(out); - - out << " MFParticleId " << GetUniqueName() << "("; - - if (fHillas) - out << "&" << fHillas->GetUniqueName(); - else - out << "\"" << fContName << "\""; - - out << ", '" << (fFilterType==kELowerThan?"<":">") << "', " << fValue << ");" << endl; - -} Index: trunk/MagicSoft/Mars/mfilter/MFAlpha.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFAlpha.h (revision 6973) +++ (revision ) @@ -1,48 +1,0 @@ -#ifndef MARS_MFAlpha -#define MARS_MFAlpha - -///////////////////////////////////////////////////////////////////////////// -// // -// MFAlpha // -// // -///////////////////////////////////////////////////////////////////////////// - -#ifndef MARS_MFilter -#include "MFilter.h" -#endif - -class MParList; -class MHillasSrc; - -class MFAlpha : public MFilter -{ -private: - MHillasSrc *fHillas; - TString fContName; - - typedef enum { kELowerThan, kEGreaterThan } FilterType_t; - FilterType_t fFilterType; - - Bool_t fResult; //! - Float_t fValue; // [deg] - - void Init(const char type, const Float_t val, - const char *name, const char *title); - - void StreamPrimitive(ofstream &out) const; - - Int_t PreProcess(MParList *pList); - Int_t Process(); - -public: - MFAlpha(const char *cname="MHillas", const char type='>', const Float_t deg=15, - const char *name=NULL, const char *title=NULL); - MFAlpha(MHillasSrc *hillas, const char type='>', const Float_t deg=15, - const char *name=NULL, const char *title=NULL); - - Bool_t IsExpressionTrue() const { return fResult; } - - ClassDef(MFAlpha, 1) // A Filter for cuts in fabs(alpha) -}; - -#endif Index: trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelBasic.cc (revision 6973) +++ (revision ) @@ -1,285 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -! Author(s): Thomas Bretz, 04/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MFSelBasic -// -// This is a class to evaluate basic cuts -// -// to be called after the calibration (when the number of photons is -// available for all pixels) -// -// The basic cuts are : -// -// remove bad runs -// thetamin < theta < thetamax -// phimin < phi < phimax -// software trigger fullfilled (with minimum no.of photons = minphotons) -// -// -///////////////////////////////////////////////////////////////////////////// - -#include "MFSelBasic.h" - -#include "MParList.h" - -#include "MPointingPos.h" - -#include "MCerPhotEvt.h" -#include "MRawRunHeader.h" - -#include "MGeomPix.h" -#include "MGeomCam.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MFSelBasic); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MFSelBasic::MFSelBasic(const char *name, const char *title) -{ - fName = name ? name : "MFSelBasic"; - fTitle = title ? title : "Filter to evaluate basic cuts"; - - // default values of cuts - SetCuts(20.0, 0.0, 60.0, 0.0, 360.0); -} - -// -------------------------------------------------------------------------- -// -// Set the cut values -// -// -void MFSelBasic::SetCuts(Float_t minphotons, - Float_t thetamin, Float_t thetamax, - Float_t phimin, Float_t phimax) -{ - fMinPhotons = minphotons; - - fThetaMin = thetamin; - fThetaMax = thetamax; - - fPhiMin = phimin; - fPhiMax = phimax; -} - -// -------------------------------------------------------------------------- -// -// Set the pointers -// -// -Int_t MFSelBasic::PreProcess(MParList *pList) -{ - fRawRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!fRawRun) - { - *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; - return kFALSE; - } - - fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); - if (!fPointPos) - { - *fLog << dbginf << "MPointingPos not found... aborting." << endl; - return kFALSE; - } - - fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); - if (!fEvt) - { - *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; - return kFALSE; - } - - fCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fCam) - { - *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; - return kFALSE; - } - - memset(fCut, 0, sizeof(fCut)); - - //------------------------- - *fLog << inf << "MFSelBasic cut values : fMinPhotons, fThetaMin, fThetaMax, fPhiMin, fPhiMax = "; - *fLog << fMinPhotons << ", " << fThetaMin << ", " << fThetaMax << ", " - << fPhiMin << ", " << fPhiMax << endl; - //------------------------- - - return kTRUE; -} - -Int_t MFSelBasic::Set(Int_t rc) -{ - fCut[rc]++; - fResult=kTRUE; - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Evaluate basic cuts -// -// bad events : fResult = kTRUE; -// good events : fResult = kFALSE; -// -Int_t MFSelBasic::Process() -{ - const Double_t theta = fPointPos->GetZd(); - const Double_t phi = fPointPos->GetAz(); - - fResult = kFALSE; - - // remove bad runs for MC gammas - //if (fMcEvt->GetEnergy() == 0.0 && fMcEvt->GetImpact() == 0.0) - //{ - // if (fRawRun->GetRunNumber() == 601 || - // fRawRun->GetRunNumber() == 613 || - // fRawRun->GetRunNumber() == 614 ) - // return Set(1); - //} - - if (thetafThetaMax) - return Set(3); - - if (phifPhiMax) - return Set(6); - - if (!SwTrigger()) - return Set(4); - - fCut[0]++; - - return kTRUE; -} -// -------------------------------------------------------------------------- -// -// Software trigger -// -// require 2 neighboring pixels (which are not in the outermost ring), -// each having at least 'fMinPhotons' photons -// -// -Bool_t MFSelBasic::SwTrigger() -{ - const Int_t entries = fEvt->GetNumPixels(); - - for (Int_t i=0; iIsPixelUsed(id1) ) - continue; - - const MGeomPix &gpix1 = (*fCam)[id1]; - if ( gpix1.IsInOutermostRing() ) - continue; - - const MCerPhotPix &pix1 = *fEvt->GetPixById(id1); - - const Double_t photons1 = pix1.GetNumPhotons(); - if (photons1 >= fMinPhotons) - return kTRUE; - } - } - return kFALSE; -} - -// -------------------------------------------------------------------------- -// -// Prints some statistics about the Basic selections. -// -Int_t MFSelBasic::PostProcess() -{ - if (GetNumExecutions()==0) - return kTRUE; - - *fLog << inf << endl; - *fLog << GetDescriptor() << " execution statistics:" << endl; - *fLog << dec << setfill(' '); - - *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ; - *fLog << (int)(fCut[1]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: bad run " << endl; - - *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) ; - *fLog << (int)(fCut[2]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: Zenith angle < " << fThetaMin << endl; - - *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) ; - *fLog << (int)(fCut[3]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: Zenith angle > " << fThetaMax << endl; - - *fLog << " " << setw(7) << fCut[5] << " (" << setw(3) ; - *fLog << (int)(fCut[5]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: Azimuth angle < " << fPhiMin << endl; - - *fLog << " " << setw(7) << fCut[6] << " (" << setw(3) ; - *fLog << (int)(fCut[6]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: Azimuth angle > " << fPhiMax << endl; - - *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) ; - *fLog << (int)(fCut[4]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: Software trigger not fullfilled" ; - *fLog << " (with fMinPhotons = " << fMinPhotons << ")" << endl; - - *fLog << " " << fCut[0] << " (" << (int)(fCut[0]*100/GetNumExecutions()) ; - *fLog << "%) Evts survived Basic selections!" << endl; - *fLog << endl; - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mfilter/MFSelBasic.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelBasic.h (revision 6973) +++ (revision ) @@ -1,73 +1,0 @@ -#ifndef MARS_MFSelBasic -#define MARS_MFSelBasic - -///////////////////////////////////////////////////////////////////////////// -// // -// MFSelBasic // -// // -// Class to evaluate basic cuts // -// // -///////////////////////////////////////////////////////////////////////////// - -#ifndef MARS_MFilter -#include "MFilter.h" -#endif - -class MPointingPos; -class MGeomCam; -class MCerPhotEvt; -//class MPedPhotCam; -class MRawRunHeader; - -class MFSelBasic : public MFilter -{ -private: - const MPointingPos *fPointPos; - const MGeomCam *fCam; // Camera Geometry - const MCerPhotEvt *fEvt; // Cerenkov Photon Event - // const MPedPhotCam *fPed; // Pedestal information - const MRawRunHeader *fRawRun; - - Float_t fMinPhotons; - Float_t fThetaMin; - Float_t fThetaMax; - - Float_t fPhiMin; - Float_t fPhiMax; - - Int_t fCut[7]; - - Bool_t fResult; - - Int_t Set(Int_t rc); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - Bool_t IsExpressionTrue() const { return fResult; } - - Bool_t SwTrigger(); - -public: - MFSelBasic(const char *name=NULL, const char *title=NULL); - - void SetCuts(Float_t minphotons=0.0, - Float_t thetamin=0.0, Float_t thetamax=90.0, - Float_t phimin=0.0, Float_t phimax=360.0); - - ClassDef(MFSelBasic, 0) // Class to evaluate basic cuts -}; - -#endif - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mfilter/MFSelFinal.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelFinal.cc (revision 6973) +++ (revision ) @@ -1,196 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MFSelFinal -// -// WHAT ARE THE FINAL CUTS? -// -// This is a class to evaluate the Final Cuts -// (these cuts define the final sample of gammas; -// relevant for the calculation of the effective collection areas) -// -// to be called after the calculation of the hadroness -// -///////////////////////////////////////////////////////////////////////////// -#include "MFSelFinal.h" - -#include - -#include "MParList.h" - -#include "MMcEvt.hxx" - -#include "MCerPhotEvt.h" - -#include "MGeomPix.h" -#include "MGeomCam.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MHillasExt.h" -#include "MHillasSrc.h" -#include "MHadronness.h" - -ClassImp(MFSelFinal); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MFSelFinal::MFSelFinal(const char *hilsrcname, - const char *name, const char *title) - : fHilSrcName(hilsrcname), fHadronnessName("MHadronness") -{ - fName = name ? name : "MFSelFinal"; - fTitle = title ? title : "Class to evaluate the Final Cuts"; - - // default values of cuts - SetCuts(1.0, 100.0, 10.); -} - -// -------------------------------------------------------------------------- -// -// Set cut values -// -// -void MFSelFinal::SetCuts(Float_t hadmax, Float_t alphamax, Float_t distmax) -{ - fHadronnessMax = hadmax; - fAlphaMax = alphamax; - fDistMax = distmax; - -} - -// -------------------------------------------------------------------------- -// -// Set the pointers -// -Int_t MFSelFinal::PreProcess(MParList *pList) -{ - fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); - if (!fHilSrc) - { - *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl; - return kFALSE; - } - - fHadronness = (MHadronness*)pList->FindObject(fHadronnessName, "MHadronness"); - if (!fHadronness) - { - *fLog << dbginf << "MHadronness not found... aborting." << endl; - return kFALSE; - } - - MGeomCam *cam = (MGeomCam*)pList->FindCreateObj("MGeomCamMagic","MGeomCam"); - if (!cam) - { - *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; - return kFALSE; - } - - fMm2Deg = cam->GetConvMm2Deg(); - - memset(fCut, 0, sizeof(fCut)); - - //---------------------- - *fLog << inf << "MFSelFinal cut values : fHadronnessMax, fAlphaMax, fDistMax = "; - *fLog << fHadronnessMax << ", " << fAlphaMax << ", " << fDistMax << endl; - //---------------------- - - return kTRUE; -} - -Int_t MFSelFinal::Set(Int_t rc) -{ - fCut[rc]++; - fResult=kTRUE; - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Evaluate final cuts -// -// if cuts are fulfilled set fResult = kTRUE -// -Int_t MFSelFinal::Process() -{ - const Double_t modalpha = fabs( fHilSrc->GetAlpha() ); - const Double_t h = fHadronness->GetHadronness(); - - fResult = kFALSE; - - if (h>fHadronnessMax) - return Set(1); - - if (modalpha>fAlphaMax) - return Set(2); - - if (fMm2Deg*fHilSrc->GetDist()>fDistMax) - return Set(3); - - fCut[0]++; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Prints some statistics about the Final selections. -// -Int_t MFSelFinal::PostProcess() -{ - if (GetNumExecutions()==0) - return kTRUE; - - *fLog << inf << endl; - *fLog << GetDescriptor() << " execution statistics:" << endl; - *fLog << dec << setfill(' '); - *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ; - *fLog << (int)(fCut[1]*100/GetNumExecutions()); - *fLog << "%) Evts skipped due to: hadronness > "<< fHadronnessMax; - *fLog << " (hadronness from '" << fHadronnessName << "')" << endl; - - *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) ; - *fLog << (int)(fCut[2]*100/GetNumExecutions()); - *fLog << "%) Evts skipped due to: |ALPHA| > " << fAlphaMax; - *fLog << " [degrees]" << endl; - - *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) ; - *fLog << (int)(fCut[3]*100/GetNumExecutions()); - *fLog << "%) Evts skipped due to: DIST > " << fDistMax; - *fLog << " [degrees]" << endl; - - *fLog << " " << fCut[0] << " (" ; - *fLog << (int)(fCut[0]*100/GetNumExecutions()); - *fLog << "%) Evts survived Final selections!" << endl; - *fLog << endl; - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mfilter/MFSelFinal.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelFinal.h (revision 6973) +++ (revision ) @@ -1,70 +1,0 @@ -#ifndef MARS_MFSelFinal -#define MARS_MFSelFinal - -///////////////////////////////////////////////////////////////////////////// -// // -// MSelFinal // -// // -// Class to evaluate final cuts // -// // -///////////////////////////////////////////////////////////////////////////// - -#ifndef MARS_MFilter -#include "MFilter.h" -#endif - -class MHillas; -class MHillasSrc; -class MHadronness; - -class MFSelFinal : public MFilter -{ -private: - MHillasSrc *fHilSrc; - MHadronness *fHadronness; - - TString fHilName; - TString fHilSrcName; - TString fHadronnessName; - - Float_t fHadronnessMax; - Float_t fAlphaMax; - Float_t fDistMax; - - Double_t fMm2Deg; // conversion mm to degrees in camera - - Int_t fCut[4]; - - Bool_t fResult; - - Int_t Set(Int_t rc); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - Bool_t IsExpressionTrue() const { return fResult; } - -public: - MFSelFinal(const char *HilSrcName="MHillasSrc", - const char *name=NULL, const char *title=NULL); - - void SetHadronnessName(const TString name) { fHadronnessName = name; } - - void SetCuts(Float_t hadmax, Float_t alphamax, Float_t distmax); - - ClassDef(MFSelFinal, 0) // Class to evaluate final cuts -}; - -#endif - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mfilter/MFSelStandard.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelStandard.cc (revision 6973) +++ (revision ) @@ -1,224 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -! Author(s): Thomas Bretz, 04/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MFSelStandard -// -// This is a class to evaluate the Standard Cuts -// -// WHAT ARE THE STANDARD CUTS? // -// -// to be called after the calculation of the image parameters -// before the g/h separation -// -///////////////////////////////////////////////////////////////////////////// - -#include "MFSelStandard.h" - -#include "MParList.h" - -#include "MMcEvt.hxx" - -#include "MGeomPix.h" -#include "MGeomCam.h" - -#include "MCerPhotEvt.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MHillas.h" -#include "MHillasExt.h" -#include "MHillasSrc.h" -#include "MNewImagePar.h" - -ClassImp(MFSelStandard); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MFSelStandard::MFSelStandard(const char *hilsrcname, - const char *name, const char *title) - : fHilName("MHillas"), fHilSrcName(hilsrcname), fImgParName("MNewImagePar") -{ - fName = name ? name : "MFSelStandard"; - fTitle = title ? title : "Class to evaluate the Standard Cuts"; - - // default values of cuts - SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0); -} - -// -------------------------------------------------------------------------- -// -// Set the values for the cuts -// -// -void MFSelStandard::SetCuts(Float_t usedpixelsmax, Float_t corepixelsmin, - Float_t sizemin, Float_t distmin, Float_t distmax, - Float_t lengthmin, Float_t widthmin) -{ - fUsedPixelsMax = usedpixelsmax; - fCorePixelsMin = corepixelsmin; - fSizeMin = sizemin; - fDistMin = distmin; - fDistMax = distmax; - fLengthMin = lengthmin; - fWidthMin = widthmin; - -} - -// -------------------------------------------------------------------------- -// -// MISSING -// -Int_t MFSelStandard::PreProcess(MParList *pList) -{ - fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); - if (!fHil) - { - *fLog << err << fHilName << " [MHillas] not found... aborting." << endl; - return kFALSE; - } - - fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); - if (!fHilSrc) - { - *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; - return kFALSE; - } - - fNewImgPar = (MNewImagePar*)pList->FindObject(fImgParName, "MNewImagePar"); - if (!fNewImgPar) - { - *fLog << err << fImgParName << " [MNewImagePar] not found... aborting." << endl; - return kFALSE; - } - - MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!cam) - { - *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; - return kFALSE; - } - - fMm2Deg = cam->GetConvMm2Deg(); - - memset(fCut, 0, sizeof(fCut)); - - //-------------------- - *fLog << inf << "MFSelStandard cut values : fUsedPixelsMax, fCorePixelsMin = "; - *fLog << fUsedPixelsMax << ", " << fCorePixelsMin << endl; - *fLog << inf << " fSizeMin, fDistMin, fDistMax = " << fSizeMin ; - *fLog << ", " << fDistMin << ", " << fDistMax << endl; - *fLog << inf << " fLengthMin, fWidthMin = " << fLengthMin ; - *fLog << ", " << fWidthMin << endl; - //-------------------- - - return kTRUE; -} - -Bool_t MFSelStandard::Set(Int_t rc) -{ - fResult = kTRUE; - fCut[rc]++; - return kTRUE; -} -// -------------------------------------------------------------------------- -// -// Evaluate standard cuts -// -// if selections are fulfilled set fResult = kTRUE; -// -// -Int_t MFSelStandard::Process() -{ - const Double_t length = fHil->GetLength() * fMm2Deg; - const Double_t width = fHil->GetWidth() * fMm2Deg; - const Double_t dist = fHilSrc->GetDist()* fMm2Deg; - //const Double_t delta = fHil->GetDelta() * kRad2Deg; - const Double_t size = fHil->GetSize(); - const Int_t numusedpixels = fNewImgPar->GetNumUsedPixels(); - const Int_t numcorepixels = fNewImgPar->GetNumCorePixels(); - - fResult = kFALSE; - - if (numusedpixels>=fUsedPixelsMax || numcorepixelsfDistMax) - return Set(3); - - if (length<=fLengthMin || width<=fWidthMin) - return Set(4); - - fCut[0]++; - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Prints some statistics about the Standard selections. -// -Int_t MFSelStandard::PostProcess() -{ - if (GetNumExecutions()==0) - return kTRUE; - - *fLog << inf << endl; - *fLog << GetDescriptor() << " execution statistics:" << endl; - *fLog << dec << setfill(' '); - *fLog << " " << setw(7) << fCut[1] << " (" << setw(3); - *fLog << (int)(fCut[1]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: Used pixels >= " << fUsedPixelsMax ; - *fLog << " or Core pixels < " << fCorePixelsMin << endl; - - *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) ; - *fLog << (int)(fCut[2]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: SIZE <= " << fSizeMin << endl; - - *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) ; - *fLog << (int)(fCut[3]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: DIST < " << fDistMin; - *fLog << " or DIST > " << fDistMax << endl; - - *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) ; - *fLog << (int)(fCut[4]*100/GetNumExecutions()) ; - *fLog << "%) Evts skipped due to: LENGTH <= " << fLengthMin; - *fLog << " or WIDTH <= " << fWidthMin << endl; - - *fLog << " " << fCut[0] << " (" ; - *fLog << (int)(fCut[0]*100/GetNumExecutions()) ; - *fLog << "%) Evts survived Standard selections!" << endl; - *fLog << endl; - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mfilter/MFSelStandard.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSelStandard.h (revision 6973) +++ (revision ) @@ -1,78 +1,0 @@ -#ifndef MARS_MFSelStandard -#define MARS_MFSelStandard - -///////////////////////////////////////////////////////////////////////////// -// // -// MFSelStandard // -// // -// Class to evaluate standard cuts // -// // -///////////////////////////////////////////////////////////////////////////// - -#ifndef MARS_MFilter -#include "MFilter.h" -#endif - -class MHillas; -class MHillasSrc; -class MNewImagePar; - -class MFSelStandard : public MFilter -{ -private: - MHillas *fHil; - MHillasSrc *fHilSrc; - MNewImagePar *fNewImgPar; - - TString fHilName; - TString fHilSrcName; - TString fImgParName; - - Float_t fUsedPixelsMax; - Float_t fCorePixelsMin; - Float_t fSizeMin; - Float_t fDistMin; - Float_t fDistMax; - Float_t fLengthMin; - Float_t fWidthMin; - - Double_t fMm2Deg; // conversion mm to degrees in camera - - Int_t fCut[5]; - - Bool_t fResult; - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - Bool_t IsExpressionTrue() const { return fResult; } - - Bool_t Set(Int_t rc); - -public: - MFSelStandard(const char *HilSrcName="MHillasSrc", - const char *name=NULL, const char *title=NULL); - - void SetHillasName(const char *name) { fHilName = name; } - void SetImgParName(const char *name) { fImgParName = name; } - - void SetCuts(Float_t usedpixelsmax, Float_t corepixelsmin, - Float_t sizemin, Float_t distmin, Float_t distmax, - Float_t lengthmin, Float_t widthmin); - - ClassDef(MFSelStandard, 0) // Class to evaluate standard cuts -}; - -#endif - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc (revision 6973) +++ (revision ) @@ -1,173 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Thomas Bretz 1/2002 -! Author(s): Wolfgang Wittek 1/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHAlphaEnergyTheta // -// // -// 3D-histogram in alpha, E-est and Theta // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHAlphaEnergyTheta.h" - -#include - -#include - -#include "MMcEvt.hxx" -#include "MHillasSrc.h" -#include "MEnergyEst.h" - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHAlphaEnergyTheta); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets name and title of the histogram. -// -MHAlphaEnergyTheta::MHAlphaEnergyTheta(const char *name, const char *title) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHAlphaEnergyTheta"; - fTitle = title ? title : "3-D histogram in alpha, energy and theta"; - - fHist.SetDirectory(NULL); - - fHist.SetTitle("3D-plot of alpha, E_{est}, Theta"); - fHist.SetXTitle("\\alpha [\\circ]"); - fHist.SetYTitle("E_{est} [GeV]"); - fHist.SetZTitle("\\Theta [\\circ]"); -} - -// -------------------------------------------------------------------------- -// -// Set binnings and prepare filling of the histogram -// -Bool_t MHAlphaEnergyTheta::SetupFill(const MParList *plist) -{ - fEnergy = (MEnergyEst*)plist->FindObject("MEnergyEst"); - if (!fEnergy) - { - *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl; - return kFALSE; - } - - fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; - return kFALSE; - } - - MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); - MBinning* binsalphaflux = (MBinning*)plist->FindObject("BinningAlphaFlux"); - MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); - if (!binsenergy || !binsalphaflux || !binstheta) - { - *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHist, binsalphaflux, binsenergy, binstheta); - - fHist.Sumw2(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the histogram -// -Bool_t MHAlphaEnergyTheta::Fill(const MParContainer *par, const Stat_t w) -{ - MHillasSrc &hil = *(MHillasSrc*)par; - - fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(), - fMcEvt->GetTelescopeTheta()*kRad2Deg, w); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHAlphaEnergyTheta::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - pad->Divide(2,2); - - TH1 *h; - - pad->cd(1); - gPad->SetBorderMode(0); - h = fHist.Project3D("expro"); - h->SetTitle("Distribution of \\alpha [\\circ]"); - h->SetXTitle("\\alpha [\\circ]"); - h->SetYTitle("Counts"); - h->Draw(opt); - h->SetBit(kCanDelete); - - pad->cd(2); - gPad->SetBorderMode(0); - gPad->SetLogx(); - h = fHist.Project3D("eypro"); - h->SetTitle("Distribution of E-est [GeV]"); - h->SetXTitle("E_{est} [GeV]"); - h->SetYTitle("Counts"); - h->Draw(opt); - h->SetBit(kCanDelete); - - pad->cd(3); - gPad->SetBorderMode(0); - h = fHist.Project3D("ezpro"); - h->SetTitle("Distribution of \\Theta [\\circ]"); - h->SetXTitle("\\Theta [\\circ]"); - h->SetYTitle("Counts"); - h->Draw(opt); - h->SetBit(kCanDelete); - - pad->cd(4); - gPad->SetBorderMode(0); - fHist.Draw(opt); - - pad->Modified(); - pad->Update(); -} Index: trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h (revision 6973) +++ (revision ) @@ -1,42 +1,0 @@ -#ifndef MARS_MHAlphaEnergyTheta -#define MARS_MHAlphaEnergyTheta - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH3 -#include -#endif - -class MMcEvt; -class MHillasSrc; -class MEnergyEst; -class MParList; - -class MHAlphaEnergyTheta : public MH -{ -private: - MMcEvt *fMcEvt; //! - MHillasSrc *fHillas; //! - MEnergyEst *fEnergy; //! - - TH3D fHist; - -public: - MHAlphaEnergyTheta(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - const TH3D *GetHist() { return &fHist; } - const TH3D *GetHist() const { return &fHist; } - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void Draw(Option_t *option=""); - - ClassDef(MHAlphaEnergyTheta, 0) //3D-histogram in alpha, Energy and theta -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.cc (revision 6973) +++ (revision ) @@ -1,279 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Thomas Bretz 1/2002 -! Author(s): Wolfgang Wittek 1/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHAlphaEnergyTime // -// // -// 3D-histogram in alpha, E-est and time // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHAlphaEnergyTime.h" - -#include - -#include - -#include "MHillasSrc.h" -#include "MEnergyEst.h" -#include "MTime.h" - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHAlphaEnergyTime); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets name and title of the histogram. -// -MHAlphaEnergyTime::MHAlphaEnergyTime(const char *name, const char *title) - : fHist() -{ - // - // set the name and title of this object - // - fName = name ? name : "MHAlphaEnergyTime"; - fTitle = title ? title : "3-D histogram in alpha, energy and time"; - - fHist.SetDirectory(NULL); - - fHist.SetTitle("3D-plot of alpha, E-est, time"); - fHist.SetXTitle("\\alpha [\\circ]"); - fHist.SetYTitle("E-est [GeV] "); - fHist.SetZTitle("time [s]"); -} - -// -------------------------------------------------------------------------- -// -// Set binnings and prepare filling of the histogram -// -Bool_t MHAlphaEnergyTime::SetupFill(const MParList *plist) -{ - fEnergy = (MEnergyEst*)plist->FindObject("MEnergyEst"); - if (!fEnergy) - { - *fLog << err << dbginf << "MHAlphaEnergyTime : MEnergyEst not found... aborting." << endl; - return kFALSE; - } - - fTime = (MTime*)plist->FindObject("MTime"); - if (!fTime) - { - *fLog << err << dbginf << "MHAlphaEnergyTime : MTime not found... aborting." << endl; - return kFALSE; - } - - MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); - MBinning* binsalphaflux = (MBinning*)plist->FindObject("BinningAlphaFlux"); - MBinning* binstime = (MBinning*)plist->FindObject("BinningTime"); - if (!binsenergy || !binsalphaflux || !binstime) - { - *fLog << err << dbginf << "MHAlphaEnergyTime : At least one MBinning not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHist, binsalphaflux, binsenergy, binstime); - - fHist.Sumw2(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the histogram -// -Bool_t MHAlphaEnergyTime::Fill(const MParContainer *par, const Stat_t w) -{ - MHillasSrc &hil = *(MHillasSrc*)par; - - fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(), *fTime, w); - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHAlphaEnergyTime::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - pad->Divide(2,2); - - TH1 *h; - - pad->cd(1); - gPad->SetBorderMode(0); - h = fHist.Project3D("ex"); - h->SetTitle("Distribution of \\alpha [\\circ]"); - h->SetXTitle("\\alpha [\\circ]"); - h->SetYTitle("Counts"); - h->Draw(opt); - h->SetBit(kCanDelete); - - pad->cd(2); - gPad->SetBorderMode(0); - gPad->SetLogx(); - h = fHist.Project3D("ey"); - h->SetTitle("Distribution of E-est [GeV]"); - h->SetXTitle("E-est [GeV] "); - h->SetYTitle("Counts"); - h->Draw(opt); - h->SetBit(kCanDelete); - - pad->cd(3); - gPad->SetBorderMode(0); - h = fHist.Project3D("ez"); - h->SetTitle("Distribution of time [s]"); - h->SetXTitle("time [s]"); - h->SetYTitle("Counts"); - h->Draw(opt); - h->SetBit(kCanDelete); - - pad->cd(4); - gPad->SetBorderMode(0); - fHist.Draw(opt); - - pad->Modified(); - pad->Update(); - -} - -// -------------------------------------------------------------------------- -// -// Integrate fHist (Alpha,E-est,Time) over the Time to get -// fAlphaEest(Alpha,E-est) -// -TH2D *MHAlphaEnergyTime::IntegrateTime(const char *title, Bool_t draw) -{ - Int_t nzbins = fHist.GetNbinsZ(); - TAxis &axez = *fHist.GetZaxis(); - axez.SetRange(1,nzbins); - - TH2D &fAlphaEest = *(TH2D *)fHist.Project3D("exy"); - - fAlphaEest.SetTitle(title); - fAlphaEest.SetXTitle("E-est [GeV] "); - fAlphaEest.SetYTitle("\\alpha [ \\circ]"); - - if (draw == kTRUE) - { - TCanvas &c = *MakeDefCanvas(title, title); - - gROOT->SetSelectedPad(NULL); - - fAlphaEest.DrawCopy(); - gPad->SetLogx(); - - c.Modified(); - c.Update(); - } - - return &fAlphaEest; -} - -// -------------------------------------------------------------------------- -// -// Integrate fHist (Alpha,E-est,Time) over E-est to get -// fAlphaTime(Alpha,Time) -// -TH2D *MHAlphaEnergyTime::IntegrateEest(const char *title, Bool_t draw) -{ - Int_t nybins = fHist.GetNbinsY(); - TAxis &axey = *fHist.GetYaxis(); - axey.SetRange(1,nybins); - - TH2D &fAlphaTime = *(TH2D *)fHist.Project3D("exz"); - - fAlphaTime.SetTitle(title); - fAlphaTime.SetXTitle("Time [s]"); - fAlphaTime.SetYTitle("\\alpha [ \\circ]"); - - if (draw == kTRUE) - { - TCanvas &c = *MakeDefCanvas(title, title); - - gROOT->SetSelectedPad(NULL); - - fAlphaTime.DrawCopy(); - - c.Modified(); - c.Update(); - } - - return &fAlphaTime; -} - -// -------------------------------------------------------------------------- -// -// Integrate fHist (Alpha,E-est,Time) over Eest and Time to get -// fAlpha(Alpha) -// -TH1D *MHAlphaEnergyTime::IntegrateEestTime(const char *title, Bool_t draw) -{ - Int_t nybins = fHist.GetNbinsY(); - TAxis &axey = *fHist.GetYaxis(); - axey.SetRange(1,nybins); - - Int_t nzbins = fHist.GetNbinsZ(); - TAxis &axez = *fHist.GetZaxis(); - axez.SetRange(1,nzbins); - - TH1D &fAlpha = *(TH1D *)fHist.Project3D("ex"); - - fAlpha.SetTitle(title); - fAlpha.SetXTitle("\\alpha [ \\circ]"); - fAlpha.SetYTitle("Counts"); - - if (draw == kTRUE) - { - TCanvas &c = *MakeDefCanvas(title, title); - - gROOT->SetSelectedPad(NULL); - - fAlpha.DrawCopy(); - - c.Modified(); - c.Update(); - } - - return &fAlpha; -} - - - - - - - Index: trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.h (revision 6973) +++ (revision ) @@ -1,57 +1,0 @@ -#ifndef MARS_MHAlphaEnergyTime -#define MARS_MHAlphaEnergyTime - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH3 -#include "TH3.h" -#endif - -#ifndef ROOT_TH2 -#include "TH2.h" -#endif - -class MHillasSrc; -class MEnergyEst; -class MTime; -class TH2D; -class MParList; - -class MHAlphaEnergyTime : public MH -{ -private: - MHillasSrc *fHillas; //! - MEnergyEst *fEnergy; //! - MTime *fTime; //! - - TH3D fHist; - -public: - MHAlphaEnergyTime(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - const TH3D *GetHist() { return &fHist; } - const TH3D *GetHist() const { return &fHist; } - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void Draw(Option_t *option=""); - - TH2D *IntegrateTime (const char *title, Bool_t Draw); - TH2D *IntegrateEest (const char *title, Bool_t Draw); - TH1D *IntegrateEestTime(const char *title, Bool_t Draw); - - ClassDef(MHAlphaEnergyTime, 0) //3D-histogram in alpha, Energy and time -}; - -#endif - - - - - - Index: trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.cc (revision 6973) +++ (revision ) @@ -1,225 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -/////////////////////////////////////////////////////////////////////// -// -// MHCT1Supercuts -// -// This class contains histograms for the supercuts -// -// the histograms are filled during the optimization of the supercuts -// -/////////////////////////////////////////////////////////////////////// -#include "MHCT1Supercuts.h" - -#include - -#include -#include -#include -#include - -#include "MParList.h" -#include "MHFindSignificance.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHCT1Supercuts); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Setup four histograms for Alpha, and Dist -// -MHCT1Supercuts::MHCT1Supercuts(const char *name, const char *title) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHCT1Supercuts"; - fTitle = title ? title : "Container for histograms for the supercuts"; - - - fDegree = new TH1F("Degree", "Degree of polynomial", 5, -0.5, 4.5); - fProb = new TH1F("Prob", "chi2 probability", 40, 0, 1.0); - fNdf = new TH1F("NDF", "NDF of polynomial fit", 60, -0.5, 59.5); - fGamma = new TH1F("gamma", "gamma", 40, 0.0, 8.0); - fNexNon = new TH1F("NexNon", "Nex / Non", 50, 0.0, 1.0); - fSigLiMa= new TH1F("Significance", "significance of gamma signal", - 60, 0.0, 120.0); - fSigtoBackg= new TH2F("SigtoBackg", "Significance vs signal/backg ratio", - 50, 0.0, 10.0, 60, 0.0, 120.0); - fSigDegree = new TH2F("SigDegree", "Significance vs Degree of polynomial", - 5, -0.5, 4.5, 60, 0.0, 120.0); - fSigNbins = new TH2F("SigNbins", "Significance vs number of bins", - 40, -0.5, 79.5, 60, 0.0, 120.0); - - fDegree->SetDirectory(NULL); - fProb->SetDirectory(NULL); - fNdf->SetDirectory(NULL); - fGamma->SetDirectory(NULL); - fNexNon->SetDirectory(NULL); - fSigLiMa->SetDirectory(NULL); - fSigtoBackg->SetDirectory(NULL); - fSigDegree->SetDirectory(NULL); - fSigNbins->SetDirectory(NULL); - - fDegree->SetXTitle("order of polynomial"); - fProb->SetXTitle("chi2 probability of polynomial fit"); - fNdf->SetXTitle("NDF of polynomial fit"); - fGamma->SetXTitle("gamma"); - fNexNon->SetXTitle("Nex / Non"); - fSigLiMa->SetXTitle("significance"); - - fSigtoBackg->SetXTitle("signa./background ratio"); - fSigtoBackg->SetYTitle("significance"); - - fSigDegree->SetXTitle("order of polynomial"); - fSigDegree->SetYTitle("significance"); - - fSigNbins->SetXTitle("number of bins"); - fSigNbins->SetYTitle("significance"); - - fDegree->SetYTitle("Counts"); - fProb->SetYTitle("Counts"); - fNdf->SetYTitle("Counts"); - fGamma->SetYTitle("Counts"); - fNexNon->SetYTitle("Counts"); - fSigLiMa->SetYTitle("Counts"); - - fSigtoBackg->SetZTitle("Counts"); - fSigDegree->SetZTitle("Counts"); - fSigNbins->SetZTitle("Counts"); -} - -// -------------------------------------------------------------------------- -// -// Delete the histograms -// -MHCT1Supercuts::~MHCT1Supercuts() -{ - delete fDegree; - delete fProb; - delete fNdf; - delete fGamma; - delete fNexNon; - delete fSigLiMa; - - delete fSigtoBackg; - delete fSigDegree; - delete fSigNbins; -} - -// -------------------------------------------------------------------------- -// -// Fill the histograms from the 'MHFindSignificance' container -// -Bool_t MHCT1Supercuts::Fill(const MParContainer *par, const Stat_t w) -{ - if (!par) - { - *fLog << err << "MHCT1Supercuts::Fill: Pointer (!=NULL) expected." << endl; - return kFALSE; - } - - MHFindSignificance &h = *(MHFindSignificance*)par; - - fDegree ->Fill(h.GetDegree( ), w); - fProb ->Fill(h.GetProb(), w); - fNdf ->Fill(h.GetNdf(), w); - fGamma ->Fill(h.GetGamma(), w); - - Double_t ratio = h.GetNon()>0.0 ? h.GetNex()/h.GetNon() : 0.0; - fNexNon ->Fill(ratio, w); - - fSigLiMa ->Fill(h.GetSigLiMa(), w); - - Double_t sigtobackg = h.GetNbg()!=0.0 ? h.GetNex() / h.GetNbg() : 0.0; - fSigtoBackg->Fill(sigtobackg, h.GetSigLiMa(), w); - - fSigDegree ->Fill(h.GetDegree(), h.GetSigLiMa(), w); - fSigNbins ->Fill(h.GetMbins(), h.GetSigLiMa(), w); - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Creates a new canvas and draws the two histograms into it. -// Be careful: The histograms belongs to this object and won't get deleted -// together with the canvas. -// -void MHCT1Supercuts::Draw(Option_t *) -{ - //TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - //pad->SetBorderMode(0); - //AppendPad(""); - - TCanvas *pad = new TCanvas("Supercuts", "Supercut plots", 900, 900); - gROOT->SetSelectedPad(NULL); - - pad->Divide(3, 3); - - pad->cd(1); - gPad->SetBorderMode(0); - fDegree->Draw(); - - pad->cd(2); - gPad->SetBorderMode(0); - fProb->Draw(); - - pad->cd(3); - gPad->SetBorderMode(0); - fNdf->Draw(); - - pad->cd(4); - gPad->SetBorderMode(0); - fGamma->Draw(); - - pad->cd(5); - gPad->SetBorderMode(0); - fNexNon->Draw(); - - pad->cd(6); - gPad->SetBorderMode(0); - fSigLiMa->Draw(); - - pad->cd(7); - gPad->SetBorderMode(0); - fSigtoBackg->Draw(); - - pad->cd(8); - gPad->SetBorderMode(0); - fSigDegree->Draw(); - - pad->cd(9); - gPad->SetBorderMode(0); - fSigNbins->Draw(); - - pad->Modified(); - pad->Update(); -} Index: trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCT1Supercuts.h (revision 6973) +++ (revision ) @@ -1,42 +1,0 @@ -#ifndef MARS_MHCT1Supercuts -#define MARS_MHCT1Supercuts - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH1F; -class TH2F; - -class MHCT1Supercuts : public MH -{ -private: - TH1F *fDegree; // order of polynomial for background fit - TH1F *fProb; // chi2 probability of polynomial fit - TH1F *fNdf; // NDF of polynomial fit - TH1F *fGamma; // Nbg = gamma * Noff - TH1F *fNexNon; // no.of excess events / no.of events in signal region - TH1F *fSigLiMa; // significance of gamma signal - - TH2F *fSigtoBackg; // significance vs signal to background ratio (Nex/Nbg) - TH2F *fSigDegree; // significance vs order of polynomial - TH2F *fSigNbins; // significance vs number of bins - - -public: - MHCT1Supercuts(const char *name=NULL, const char *title=NULL); - ~MHCT1Supercuts(); - - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - void Draw(Option_t *opt=NULL); - - ClassDef(MHCT1Supercuts, 1) // Container which holds histograms for the supercuts -}; - -#endif - - - - - Index: trunk/MagicSoft/Mars/mhist/MHCurrents.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCurrents.cc (revision 6973) +++ (revision ) @@ -1,246 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 12/2002 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MHCurrents -// -///////////////////////////////////////////////////////////////////////////// -#include "MHCurrents.h" - -#include - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" -#include "MBinning.h" -#include "MCurrents.h" -#include "MCamDisplay.h" - -#include "MGeomCam.h" -#include "MGeomPix.h" - -ClassImp(MHCurrents); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Reset all pixels to 0 and reset fEntries to 0. -// -void MHCurrents::Clear(const Option_t *) -{ - const Int_t n = fCam ? fCam->GetNumPixels() : 577; - - // FIXME: Implement a clear function with setmem - fSum.Set(n); // also clears memory - fRms.Set(n); -/* for (int i=0; i<577; i++) - { - fSum[i] = 0; - fRms[i] = 0; - }*/ - - fEntries = 0; -} - -// -------------------------------------------------------------------------- -// -// Initialize the name and title of the task. -// Resets the sum histogram -// -MHCurrents::MHCurrents(const char *name, const char *title) - : /*fSum(577), fRms(577), */fCam(NULL), fEvt(NULL), fDispl(NULL) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHCurrents"; - fTitle = title ? title : "Average of MCurrents"; - - Clear(); - - fHist.SetName("currents;avg"); - fHist.SetTitle("Avg.Currents [nA]"); - fHist.SetXTitle("Pixel Index"); - fHist.SetYTitle("A [nA]"); - fHist.SetDirectory(NULL); - fHist.SetLineColor(kGreen); - fHist.SetMarkerStyle(kFullDotMedium); - fHist.SetMarkerSize(0.3); -} - -// -------------------------------------------------------------------------- -// -// Delete the corresponding camera display if available -// -MHCurrents::~MHCurrents() -{ - if (fDispl) - delete fDispl; -} - -// -------------------------------------------------------------------------- -// -// Get the event (MCerPhotEvt) the histogram might be filled with. If -// it is not given, it is assumed, that it is filled with the argument -// of the Fill function. -// Looks for the camera geometry MGeomCam and resets the sum histogram. -// -Bool_t MHCurrents::SetupFill(const MParList *plist) -{ - fEvt = (MCurrents*)plist->FindObject("MCurrents"); - if (!fEvt) - *fLog << warn << GetDescriptor() << ": No MCerPhotEvt available..." << endl; - - fCam = (MGeomCam*)plist->FindObject("MGeomCam"); - /* - if (!fCam) - *fLog << warn << GetDescriptor() << ": No MGeomCam found... assuming Magic geometry!" << endl; - */ - if (!fCam) - { - *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl; - return kFALSE; - } - - Clear(); - - const Int_t n = fSum.GetSize(); - - MBinning bins; - bins.SetEdges(n, -0.5, n-0.5); - bins.Apply(fHist); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the histograms with data from a MCerPhotEvt-Container. -// -Bool_t MHCurrents::Fill(const MParContainer *par, const Stat_t w) -{ - const MCurrents *evt = par ? (MCurrents*)par : fEvt; - if (!evt) - { - *fLog << err << dbginf << "No MCurrents found..." << endl; - return kFALSE; - } - - const Int_t n = fSum.GetSize(); - for (UInt_t idx=0; idxGetPixelContent(val, idx)) - continue; - - fSum[idx] += val; - fRms[idx] += val*val; - - } - - fEntries++; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Scale the sum container with the number of entries -// -Bool_t MHCurrents::Finalize() -{ - if (fEntries<2) - { - *fLog << warn << "WARNING - " << GetDescriptor() << " doesn't contain enough entries." << endl; - return kTRUE; - } - - const Int_t n = fSum.GetSize(); - for (UInt_t i=0; i was set to 0 " << endl; - fRms[i]=0; - } - - else - fRms[i] = TMath::Sqrt(fRms[i]); - - // calc mean value for pixel index i - fSum[i] /= fEntries; - - fHist.SetBinContent(i+1, fSum[i]); - fHist.SetBinError( i+1, fRms[i]); - } - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the present 'fill status' -// -void MHCurrents::Draw(Option_t *o) -{ - if (!fCam) - { - *fLog << warn << "WARNING - Cannot draw " << GetDescriptor() << ": No Camera Geometry available." << endl; - return; - } - - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this, 750, 600); - pad->SetBorderMode(0); - - SetDrawOption(o); - AppendPad(""); -} - -// -------------------------------------------------------------------------- -// -// If a camera display is not yet assigned, assign a new one. -// -void MHCurrents::Paint(Option_t *option) -{ - if (!fCam) - { - *fLog << warn << "WARNING - Cannot paint " << GetDescriptor() << ": No Camera Geometry available." << endl; - return; - } - - if (!fDispl) - fDispl = new MCamDisplay(fCam); - - TString opt(GetDrawOption()); - - fDispl->Fill(opt.Contains("rms", TString::kIgnoreCase) ? fRms : fSum); - fDispl->Paint(); -} Index: trunk/MagicSoft/Mars/mhist/MHCurrents.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCurrents.h (revision 6973) +++ (revision ) @@ -1,55 +1,0 @@ -#ifndef MARS_MHCurrents -#define MARS_MHCurrents - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH1 -#include -#endif - -#ifndef ROOT_TArrayF -#include -#endif - -class MCurrents; -class MGeomCam; -class MCamDisplay; - -class MHCurrents : public MH -{ -private: - TArrayF fSum; // storing the sum - TArrayF fRms; // storing the rms - Int_t fEntries; // number of entries in the histogram - MGeomCam *fCam; // the present geometry - MCurrents *fEvt; //! the current event - MCamDisplay *fDispl; //! the camera display - - TH1F fHist; - -public: - MHCurrents(const char *name=NULL, const char *title=NULL); - ~MHCurrents(); - - void Clear(const Option_t *o=""); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - Bool_t Finalize(); - - TH1 *GetHistByName(const TString name) { return NULL; } - - const TArrayF &GetSum() const { return fSum; } - const TArrayF &GetRms() const { return fRms; } - - const TH1F &GetHist() const { return fHist; } - - void Draw(Option_t *opt=""); - void Paint(Option_t *option=""); - - ClassDef(MHCurrents, 1) // Histogram to sum camera events -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/MHEnergyTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEnergyTheta.cc (revision 6973) +++ (revision ) @@ -1,133 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Thomas Bretz 1/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHEnergyTheta // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHEnergyTheta.h" - -#include - -#include "MH.h" -#include "MBinning.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MMcEvt.hxx" - -#include "MParList.h" - -ClassImp(MHEnergyTheta); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Creates the three necessary histograms: -// - selected showers (input) -// - all showers (input) -// - collection area (result) -// -MHEnergyTheta::MHEnergyTheta(const char *name, const char *title) -{ - // initialize the histogram for the distribution r vs E - // - // we set the energy range from 1 Gev to 10000 GeV (in log 5 orders - // of magnitude) and for each order we take 10 subdivision --> 50 xbins - // - // we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins - - - fName = name ? name : "MHEnergyTheta"; - fTitle = title ? title : "Data to Calculate Collection Area"; - - fHist.SetDirectory(NULL); - - fHist.SetXTitle("E [GeV]"); - fHist.SetYTitle("\\Theta [\\circ]"); - fHist.SetZTitle("N"); -} - -Bool_t MHEnergyTheta::SetupFill(const MParList *plist) -{ - MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); - MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); - if (!binsenergy || !binstheta) - { - *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHist, binsenergy, binstheta); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill data into the histogram which contains all showers -// -Bool_t MHEnergyTheta::Fill(const MParContainer *par, const Stat_t w) -{ - const MMcEvt &mcevt = *(MMcEvt*)par; - - fHist.Fill(mcevt.GetEnergy(), mcevt.GetTheta(), w); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram with all showers -// -void MHEnergyTheta::Draw(Option_t* option) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - pad->SetLogy(); - - fHist.Draw(option); - - pad->Modified(); - pad->Update(); -} - -// -------------------------------------------------------------------------- -// -// Calculate the Efficiency (collection area) and set the 'ReadyToSave' -// flag -// -void MHEnergyTheta::Divide(const TH2D *h1, const TH2D *h2) -{ - // Description! - - fHist.Sumw2(); - fHist.Divide((TH2D*)h1, (TH2D*)h2); - - SetReadyToSave(); -} Index: trunk/MagicSoft/Mars/mhist/MHEnergyTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEnergyTheta.h (revision 6973) +++ (revision ) @@ -1,42 +1,0 @@ -#ifndef MARS_MHEnergyTheta -#define MARS_MHEnergyTheta - -#ifndef MARS_MH -#include "MH.h" -#endif -#ifndef ROOT_TH2 -#include -#endif - -class MTime; -class MParList; - -class MHEnergyTheta : public MH -{ -private: - TH2D fHist; //! - -public: - MHEnergyTheta(const char *name=NULL, const char *title=NULL); - - Bool_t Fill(const MParContainer *cont, const Stat_t w=1); - - void Draw(Option_t *option=""); - - Bool_t SetupFill(const MParList *plist); - - const TH2D *GetHist() { return &fHist; } - const TH2D *GetHist() const { return &fHist; } - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void Divide(const TH2D *h1, const TH2D *h2); - void Divide(const MHEnergyTheta *h1, const MHEnergyTheta *h2) - { - Divide(h1->GetHist(), h2->GetHist()); - } - - ClassDef(MHEnergyTheta, 1) // Histogram: Energy vs. Theta -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/MHEnergyTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEnergyTime.cc (revision 6973) +++ (revision ) @@ -1,142 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Thomas Bretz 1/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHEnergyTime // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHEnergyTime.h" - -#include - -#include "MMcEvt.hxx" -#include "MTime.h" - -#include "MH.h" -#include "MBinning.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" - -ClassImp(MHEnergyTime); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Creates the three necessary histograms: -// - selected showers (input) -// - all showers (input) -// - collection area (result) -// -MHEnergyTime::MHEnergyTime(const char *name, const char *title) -{ - // initialize the histogram for the distribution r vs E - // - // we set the energy range from 1 Gev to 10000 GeV (in log 5 orders - // of magnitude) and for each order we take 10 subdivision --> 50 xbins - // - // we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins - - fName = name ? name : "MHEnergyTime"; - fTitle = title ? title : "Data to Calculate Collection Area"; - - fHist.SetDirectory(NULL); - - fHist.SetXTitle("E [GeV]"); - fHist.SetYTitle("t [s]"); - fHist.SetZTitle("N"); -} - -Bool_t MHEnergyTime::SetupFill(const MParList *plist) -{ - fTime = (MTime*)plist->FindObject("MTime"); - if (!fTime) - { - *fLog << err << dbginf << "MTime not found... aborting." << endl; - return kFALSE; - } - - const MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); - const MBinning* binstime = (MBinning*)plist->FindObject("BinningTime"); - if (!binsenergy || !binstime) - { - *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHist, binsenergy, binstime); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill data into the histogram which contains all showers -// -Bool_t MHEnergyTime::Fill(const MParContainer *par, const Stat_t w) -{ - const MMcEvt &mcevt = *(MMcEvt*)par; - - fHist.Fill(mcevt.GetEnergy(), *fTime, w); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram with all showers -// -void MHEnergyTime::Draw(Option_t* option) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - pad->SetLogy(); - fHist.DrawCopy(option); - - pad->Modified(); - pad->Update(); -} - -// -------------------------------------------------------------------------- -// -// Calculate the Efficiency (collection area) and set the 'ReadyToSave' -// flag -// -void MHEnergyTime::Divide(const TH2D *h1, const TH2D *h2) -{ - // Description! - - fHist.Sumw2(); - fHist.Divide((TH2D*)h1, (TH2D*)h2); - - SetReadyToSave(); -} Index: trunk/MagicSoft/Mars/mhist/MHEnergyTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEnergyTime.h (revision 6973) +++ (revision ) @@ -1,43 +1,0 @@ -#ifndef MARS_MHEnergyTime -#define MARS_MHEnergyTime - -#ifndef MARS_MH -#include "MH.h" -#endif -#ifndef ROOT_TH2 -#include "TH2.h" -#endif - -class MTime; -class MParList; - -class MHEnergyTime : public MH -{ -private: - MTime *fTime; - TH2D fHist; //! - -public: - MHEnergyTime(const char *name=NULL, const char *title=NULL); - - Bool_t Fill(const MParContainer *cont, const Stat_t w=1); - - void Draw(Option_t *option=""); - - Bool_t SetupFill(const MParList *plist); - - const TH2D *GetHist() { return &fHist; } - const TH2D *GetHist() const { return &fHist; } - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void Divide(const TH2D *h1, const TH2D *h2); - void Divide(const MHEnergyTime *h1, const MHEnergyTime *h2) - { - Divide(h1->GetHist(), h2->GetHist()); - } - - ClassDef(MHEnergyTime, 1) // Histogram: Energy vs. Time -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/MHSigmaPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaPixel.cc (revision 6973) +++ (revision ) @@ -1,138 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Robert Wagner 10/2002 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHSigmaPixel // -// // -// 2D-Histogram pedestal sigma vs pixel number // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHSigmaPixel.h" - -#include - -#include - -#include "MPedestalCam.h" -#include "MPedestalPix.h" -#include "MSigmabar.h" - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHSigmaPixel); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets name and title of the histogram. -// -MHSigmaPixel::MHSigmaPixel(const char *name, const char *title) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHSigmaPixel"; - fTitle = title ? title : "2-D histogram in sigma and pixel"; - - fHist.SetDirectory(NULL); - - fHist.SetTitle("Sigma vs pixel #"); - fHist.SetXTitle("pixel #"); - fHist.SetYTitle("\\sigma"); - fHist.SetZTitle("N"); -} - -// -------------------------------------------------------------------------- -// -// Set binnings and prepare filling of the histogram. The binning for the -// pixel axis is derived automagically from the MPedestalCam container -// expected to be found in the MParList *plist. -// -Bool_t MHSigmaPixel::SetupFill(const MParList *plist) -{ - fPedestalCam = (MPedestalCam*)plist->FindObject("MPedestalCam"); - if (!fPedestalCam) - { - *fLog << err << "MPedestalCam not found... aborting." << endl; - return kFALSE; - } - - MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigma"); - if (!binssigma) - { - *fLog << err << "BinningSigma [MBinning] not found... aborting." << endl; - return kFALSE; - } - - const Int_t n = fPedestalCam->GetSize(); - - MBinning binspixel; - binspixel.SetEdges(n, -0.5, -0.5+n); - - SetBinning(&fHist, &binspixel, binssigma); - - fHist.Sumw2(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the histogram -// -Bool_t MHSigmaPixel::Fill(const MParContainer *par, const Stat_t w) -{ - const MPedestalCam &ped = *(MPedestalCam*)par; - for (Int_t i=0;i<(ped.GetSize());i++) - { - const MPedestalPix pix = ped[i]; - fHist.Fill(i, pix.GetPedestalRms()); - } - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHSigmaPixel::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - fHist.Draw(opt); - - gPad->Modified(); - gPad->Update(); -} - Index: trunk/MagicSoft/Mars/mhist/MHSigmaPixel.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaPixel.h (revision 6973) +++ (revision ) @@ -1,46 +1,0 @@ -#ifndef MARS_MHSigmaPixel -#define MARS_MHSigmaPixel - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH2 -#include -#endif - -class TH2D; - -class MMcEvt; -class MParList; -class MPedestalCam; - -class MHSigmaPixel : public MH -{ -private: - MPedestalCam *fPedestalCam; //! - - TH2D fHist; - -public: - MHSigmaPixel(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - const TH2D *GetHist() { return &fHist; } - const TH2D *GetHist() const { return &fHist; } - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void Draw(Option_t *option=""); - - ClassDef(MHSigmaPixel, 1) //2D-histogram in Sigma and Pixel number -}; - -#endif - - - - - Index: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 6973) +++ (revision ) @@ -1,512 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 1/2003 -! Author(s): Thomas Bretz 4/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHSigmaTheta // -// // -// calculates - the 2D-histogram sigmabar vs. Theta (Inner), and // -// - the 2D-histogram sigmabar vs. Theta (Outer) // -// - the 3D-histogram sigma, pixel no., Theta // -// - the 3D-histogram (sigma^2-sigmabar^2), pixel no., Theta // -// - the 2D-histogram Theta vs. Phi // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHSigmaTheta.h" - -#include - -#include "MTime.h" -#include "MPointingPos.h" - -#include "MBinning.h" -#include "MParList.h" - -#include "MGeomCam.h" -#include "MGeomPix.h" -#include "MBadPixelsCam.h" - -#include "MPedPhotCam.h" -#include "MPedPhotPix.h" - -#include "MCerPhotEvt.h" -#include "MCerPhotPix.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHSigmaTheta); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Constructor. -// -MHSigmaTheta::MHSigmaTheta(const char *name, const char *title) -{ - fName = name ? name : "MHSigmaTheta"; - fTitle = title ? title : "2D histogram sigmabar vs. Theta"; - - fThetaPhi = new TH2D; - fThetaPhi->SetDirectory(NULL); - fThetaPhi->UseCurrentStyle(); - fThetaPhi->SetName("2D-ThetaPhi"); - fThetaPhi->SetTitle("2D: \\Theta vs. \\Phi"); - fThetaPhi->SetXTitle("\\phi [\\circ]"); - fThetaPhi->SetYTitle("\\Theta [\\circ]"); - fThetaPhi->SetTitleOffset(1.2,"Y"); - - fSigmaTheta = new TH2D; - fSigmaTheta->SetDirectory(NULL); - fSigmaTheta->UseCurrentStyle(); - fSigmaTheta->SetName("2D-ThetaSigmabar(Inner)"); - fSigmaTheta->SetTitle("2D: \\bar{\\sigma}, \\Theta"); - fSigmaTheta->SetXTitle("\\Theta [\\circ]"); - fSigmaTheta->SetYTitle("Sigmabar(Inner)"); - fSigmaTheta->SetTitleOffset(1.2,"Y"); - - fSigmaThetaOuter = new TH2D; - fSigmaThetaOuter->SetDirectory(NULL); - fSigmaThetaOuter->UseCurrentStyle(); - fSigmaThetaOuter->SetName("2D-ThetaSigmabar(Outer)"); - fSigmaThetaOuter->SetTitle("2D: \\bar{\\sigma}, \\Theta"); - fSigmaThetaOuter->SetXTitle("\\Theta [\\circ]"); - fSigmaThetaOuter->SetYTitle("Sigmabar(Outer)"); - fSigmaThetaOuter->SetTitleOffset(1.2,"Y"); - - fSigmaPixTheta = new TH3D; - fSigmaPixTheta->SetDirectory(NULL); - fSigmaPixTheta->UseCurrentStyle(); - fSigmaPixTheta->SetName("3D-ThetaPixSigma"); - fSigmaPixTheta->SetTitle("3D: \\Theta, Pixel Id, \\sigma"); - fSigmaPixTheta->SetXTitle("\\Theta [\\circ]"); - fSigmaPixTheta->SetYTitle("Pixel Id"); - fSigmaPixTheta->SetZTitle("Sigma"); - - fDiffPixTheta = new TH3D; - fDiffPixTheta->SetDirectory(NULL); - fDiffPixTheta->UseCurrentStyle(); - fDiffPixTheta->SetName("3D-ThetaPixDiff"); - //fDiffPixTheta->SetTitle("3D: \\Theta, Pixel Id, {\\sigma}^{2}-\\bar{\\sigma}^{2}"); - fDiffPixTheta->SetTitle("3D: \\Theta, Pixel Id, (Sigma2-Sigmabar2)/Area"); - fDiffPixTheta->SetXTitle("\\Theta [\\circ]"); - fDiffPixTheta->SetYTitle("Pixel Id"); - fDiffPixTheta->SetZTitle("(Sigma2 - Sigmabar2)/Area"); - - // Define default binning - fBinsPhi = new MBinning; - fBinsPhi->SetEdges( 20, 0.0, 360.0); - - Double_t fThetaLo = 0.0; - Double_t fThetaHi = 90.0; - fBinsTheta = new MBinning; - fBinsTheta->SetEdgesCos( 10, fThetaLo, fThetaHi); // theta - //fBinsTheta->SetEdges( 1, fThetaLo, fThetaHi); // theta - - fBinsSigma = new MBinning; - fBinsSigma->SetEdges( 100, 0.0, 120.0); // sigma - - fBinsSigmabarIn = new MBinning; - fBinsSigmabarOut = new MBinning;; - fBinsSigmabarIn->SetEdges( 100, 0.0, 25.0); // sigmabar (inner) - fBinsSigmabarOut->SetEdges(100, 0.0, 60.0); // sigmabar (outer) - - fBinsPix = new MBinning; - fBinsPix->SetEdges(578, -0.5, 577.5); - - fBinsDiff = new MBinning; - fBinsDiff->SetEdges( 100, -500.0, 1500.0); // (sigma2-sigmabar2)/area - - SetBinning(fThetaPhi, fBinsPhi, fBinsTheta); - SetBinning(fSigmaTheta, fBinsTheta, fBinsSigmabarIn); - SetBinning(fSigmaThetaOuter, fBinsTheta, fBinsSigmabarOut); - SetBinning(fSigmaPixTheta, fBinsTheta, fBinsPix, fBinsSigma); - SetBinning(fDiffPixTheta, fBinsTheta, fBinsPix, fBinsDiff); - - //-------------------------------------------- - fNamePedPhotCam = "MPedPhotCamFromData"; -} - - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MHSigmaTheta::~MHSigmaTheta() -{ - delete fThetaPhi; - delete fSigmaTheta; - delete fSigmaThetaOuter; - delete fSigmaPixTheta; - delete fDiffPixTheta; - - delete fBinsPhi; - delete fBinsTheta; - delete fBinsSigma; - delete fBinsSigmabarIn; - delete fBinsSigmabarOut; - delete fBinsPix; - delete fBinsDiff; -} - -// -------------------------------------------------------------------------- -// -// Set the binnings and prepare the filling of the histogram -// -Bool_t MHSigmaTheta::SetupFill(const MParList *plist) -{ - fCam = (MGeomCam*)plist->FindObject("MGeomCam"); - if (!fCam) - { - *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl; - return kFALSE; - } - - fPointPos = (MPointingPos*)plist->FindObject("MPointingPos"); - if (!fPointPos) - { - *fLog << err << "MPointingPos not found... aborting." << endl; - return kFALSE; - } - - fPed = (MPedPhotCam*)plist->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam"); - if (!fPed) - { - *fLog << err << AddSerialNumber(fNamePedPhotCam) - << "[MPedPhotCam] not found... aborting." << endl; - return kFALSE; - } - fPed->InitSize(fCam->GetNumPixels()); - - - fBad = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam"); - if (!fBad) - { - *fLog << err << "MBadPixelsCam not found... continue. " << endl; - } - - - fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt"); - if (!fEvt) - { - *fLog << err << "MCerPhotEvt not found... aborting." << endl; - return kFALSE; - } - - - //--------------------------------------------------------------- - *fLog << inf << "Name of MPedPhotCam container : " - << fNamePedPhotCam << endl; - - - //--------------------------------------------------------------- - - // Get Phi Binning - MBinning* binsphi = (MBinning*)plist->FindObject("BinningPhi", "MBinning"); - if (!binsphi) - { - *fLog << warn << "Object 'BinningPhi' [MBinning] not found... use default binning." << endl; - binsphi = fBinsPhi; - } - - // Get Theta Binning - MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning"); - if (!binstheta) - { - *fLog << warn << "Object 'BinningTheta' [MBinning] not found... use default binning." << endl; - binstheta = fBinsTheta; - } - - // Get Sigma binning - MBinning* binssig = (MBinning*)plist->FindObject("BinningSigma", "MBinning"); - if (!binssig) - { - *fLog << warn << "Object 'BinningSigma' [MBinning] not found... use default binning." << endl; - binssig = fBinsSigma; - } - - // Get SigmabarIn binning - MBinning* binssigmain = (MBinning*)plist->FindObject("BinningSigmabarIn", "MBinning"); - if (!binssigmain) - { - *fLog << warn << "Object 'BinningSigmabarIn' [MBinning] not found... use default binning." << endl; - binssigmain = fBinsSigmabarIn; - } - - // Get SigmabarOut binning - MBinning* binssigmaout = (MBinning*)plist->FindObject("BinningSigmabarOut", "MBinning"); - if (!binssigmaout) - { - *fLog << warn << "Object 'BinningSigmabarOut' [MBinning] not found... use default binning." << endl; - binssigmaout = fBinsSigmabarOut; - } - - // Get binning for (sigma^2-sigmabar^2) - MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2", "MBinning"); - if (!binsdiff) - { - *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... use default binning." << endl; - binsdiff = fBinsDiff; - } - - - //--------------------------------------------------------------- - - // Get binning for pixel number - const UInt_t npix1 = fPed->GetSize()+1; - //*fLog << "MHSigmaTheata::SetupFill(); npix1 = " << npix1 << endl; - MBinning binspix("BinningPixel"); - binspix.SetEdges(npix1, -0.5, npix1-0.5); - - // Set binnings in histograms - SetBinning(fThetaPhi, binsphi, binstheta); - SetBinning(fSigmaTheta, binstheta, binssigmain); - SetBinning(fSigmaThetaOuter, binstheta, binssigmaout); - SetBinning(fSigmaPixTheta, binstheta, &binspix, binssig); - SetBinning(fDiffPixTheta, binstheta, &binspix, binsdiff); - - *fLog << "MHSigmaTheta::SetupFill(); binnings were set" << endl; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the histograms -// -// ignore pixels if they are unused or blind -// -Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w) -{ - Double_t theta = fPointPos->GetZd(); - Double_t phi = fPointPos->GetAz(); - - fPed->ReCalc(*fCam, fBad); - - Double_t mysig = (fPed->GetArea(0)).GetRms(); - Double_t mysigouter = (fPed->GetArea(1)).GetRms(); - - //*fLog << "theta, mysig, mysigouter = " << theta << ", " << mysig - // << ", " << mysigouter << endl; - - fThetaPhi->Fill(phi, theta); - fSigmaTheta->Fill(theta, mysig); - fSigmaThetaOuter->Fill(theta, mysigouter); - - MCerPhotPix *cerpix = 0; - TIter Next(*fEvt); - while ( (cerpix=(MCerPhotPix*)Next()) ) - { - const Int_t id = cerpix->GetPixId(); - - if (!cerpix->IsPixelUsed()) - { - //*fLog << all << "MHSigmaTheta::Fill; unused pixel found, id = " - // << id << endl; - continue; - } - - const MPedPhotPix &pix = (*fPed)[id]; - - // ratio is the area of pixel 0 - // divided by the area of the current pixel - const Double_t ratio = fCam->GetPixRatio(id); - const Double_t sigma = pix.GetRms(); - - fSigmaPixTheta->Fill(theta, (Double_t)id, sigma); - - Double_t diff; - const Byte_t aidx = (*fCam)[id].GetAidx(); - if (aidx == 0) - { - // inner pixel - diff = (sigma*sigma - mysig*mysig) * ratio; - } - else - { - // outer pixel - diff = (sigma*sigma - mysigouter*mysigouter) * ratio; - } - - fDiffPixTheta->Fill(theta, (Double_t)id, diff); - } - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Update the projections and (if possible) set log scales before painting -// -void MHSigmaTheta::Paint(Option_t *opt) -{ - TVirtualPad *padsave = gPad; - - TH1D* h; - - padsave->cd(1); - if ((h = (TH1D*)gPad->FindObject("ProjX-Theta"))) - { - ProjectionX(*h, *fSigmaTheta); - if (h->GetEntries()!=0) - gPad->SetLogy(); - } - - padsave->cd(4); - if ((h = (TH1D*)gPad->FindObject("ProjY-sigma"))) - ProjectionY(*h, *fSigmaTheta); - - gPad = padsave; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHSigmaTheta::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - AppendPad(""); - - pad->Divide(3, 3); - - // draw the 2D histogram Sigmabar versus Theta - TH1 *h; - - pad->cd(1); - gPad->SetBorderMode(0); - h = fSigmaTheta->ProjectionX("ProjX-Theta", -1, 9999, "E"); - h->SetDirectory(NULL); - h->UseCurrentStyle(); - h->SetTitle("Distribution of \\Theta"); - h->SetXTitle("\\Theta [\\circ]"); - h->SetYTitle("No.of events"); - h->SetTitleOffset(1.2,"Y"); - h->Draw(opt); - h->SetBit(kCanDelete); - - pad->cd(2); - gPad->SetBorderMode(0); - h = fDiffPixTheta->Project3D("zx"); - h->SetDirectory(NULL); - h->UseCurrentStyle(); - //h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)"); - h->SetTitle("(Sigma2-Sigmabar2)/Area vs. \\Theta (all pixels)"); - h->SetXTitle("\\Theta [\\circ]"); - h->SetYTitle("(Sigma2 - Sigmabar2) / Area"); - h->SetTitleOffset(1.2,"Y"); - h->Draw("box"); - h->SetBit(kCanDelete); - - pad->cd(3); - gPad->SetBorderMode(0); - h = fSigmaPixTheta->Project3D("zx"); - h->SetDirectory(NULL); - h->UseCurrentStyle(); - h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)"); - h->SetXTitle("\\Theta [\\circ]"); - h->SetYTitle("Sigma"); - h->SetTitleOffset(1.2,"Y"); - h->Draw("box"); - h->SetBit(kCanDelete); - - //pad->cd(7); - //gPad->SetBorderMode(0); - //h = fSigmaTheta->ProjectionY("ProjY-sigma", -1, 9999, "E"); - //h->SetDirectory(NULL); - //h->UseCurrentStyle(); - //h->SetTitle("Distribution of \\bar{\\sigma}_{ped}"); - //h->SetXTitle("\\bar{\\sigma}_{ped}"); - //h->SetYTitle("No.of events"); - //h->SetTitleOffset(1.2,"Y"); - //h->Draw(opt); - //h->SetBit(kCanDelete); - - pad->cd(5); - gPad->SetBorderMode(0); - h = fDiffPixTheta->Project3D("zy"); - h->SetDirectory(NULL); - h->UseCurrentStyle(); - //h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)"); - h->SetTitle("(Sigma2-Sigmabar2)/Area vs. pixel Id (all \\Theta)"); - h->SetXTitle("Pixel Id"); - h->SetYTitle("(Sigma2 - sigmabar2) / Area"); - h->SetTitleOffset(1.2,"Y"); - h->Draw("box"); - h->SetBit(kCanDelete); - - pad->cd(6); - gPad->SetBorderMode(0); - h = fSigmaPixTheta->Project3D("zy"); - h->SetDirectory(NULL); - h->UseCurrentStyle(); - h->SetTitle("\\sigma_{ped} vs. pixel Id (all \\Theta)"); - h->SetXTitle("Pixel Id"); - h->SetYTitle("Sigma"); - h->SetTitleOffset(1.2,"Y"); - h->Draw("box"); - h->SetBit(kCanDelete); - - pad->cd(4); - gPad->SetBorderMode(0); - fSigmaTheta->Draw(); - - pad->cd(7); - gPad->SetBorderMode(0); - fSigmaThetaOuter->Draw(); - - pad->cd(8); - gPad->SetBorderMode(0); - fThetaPhi->Draw(); - - pad->cd(9); - gPad->SetBorderMode(0); - h = fThetaPhi->ProjectionX("ProjX-Phi", -1, 9999, "E"); - h->SetDirectory(NULL); - h->UseCurrentStyle(); - h->SetTitle("Distribution of \\Phi"); - h->SetXTitle("\\phi [\\circ]"); - h->SetYTitle("No. of events"); - h->SetTitleOffset(1.2,"Y"); - h->Draw(opt); - h->SetBit(kCanDelete); -} - - - - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 6973) +++ (revision ) @@ -1,89 +1,0 @@ -#ifndef MARS_MHSigmaTheta -#define MARS_MHSigmaTheta - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH2 -#include -#endif - -#ifndef ROOT_TH3 -#include -#endif - - -class MParList; -class MGeomCam; -class MCerPhotEvt; -class MPointingPos; -class MPedPhotCam; -class MBadPixelsCam; - - -class MHSigmaTheta : public MH -{ -private: - const MGeomCam *fCam; //! - MPedPhotCam *fPed; //! - MCerPhotEvt *fEvt; //! - MPointingPos *fPointPos; //! - MBadPixelsCam *fBad; //! - - TString fNamePedPhotCam; //! name of the 'MPedPhotCam' container - - // sigmabar is the average pedestal sigma - TH2D *fSigmaTheta; // 2D-distribution sigmabar versus Theta (Inner) - TH2D *fSigmaThetaOuter; // 2D-distribution sigmabar versus Theta (Outer) - - TH3D *fSigmaPixTheta; // 3D-distr.:Theta, pixel, pedestal sigma - TH3D *fDiffPixTheta; // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2 - - TH2D *fThetaPhi; // 2D-distribution Theta versus Phi - - MBinning *fBinsPhi; //! - MBinning *fBinsTheta; //! - MBinning *fBinsSigma; //! - MBinning *fBinsSigmabarIn; //! - MBinning *fBinsSigmabarOut; //! - MBinning *fBinsPix; //! - MBinning *fBinsDiff; //! - - void Paint(Option_t *opt=""); - -public: - MHSigmaTheta(const char *name=NULL, const char *title=NULL); - ~MHSigmaTheta(); - - void SetNamePedPhotCam(const char *name) { fNamePedPhotCam = name; } - - Bool_t SetupFill(const MParList *plist); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - const TH2D *GetThetaPhi() { return fThetaPhi; } - const TH2D *GetThetaPhi() const { return fThetaPhi; } - - const TH2D *GetSigmaTheta() { return fSigmaTheta; } - const TH2D *GetSigmaTheta() const { return fSigmaTheta; } - - const TH2D *GetSigmaThetaOuter() { return fSigmaThetaOuter; } - const TH2D *GetSigmaThetaOuter() const { return fSigmaThetaOuter; } - - const TH3D *GetSigmaPixTheta() { return fSigmaPixTheta; } - const TH3D *GetSigmaPixTheta() const { return fSigmaPixTheta; } - - const TH3D *GetDiffPixTheta() { return fDiffPixTheta; } - const TH3D *GetDiffPixTheta() const { return fDiffPixTheta; } - - void Draw(Option_t *option=""); - - ClassDef(MHSigmaTheta, 1) //2D-histogram sigmabar vs. Theta -}; - -#endif - - - - - Index: trunk/MagicSoft/Mars/mhist/MHSigmabarTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmabarTheta.cc (revision 6973) +++ (revision ) @@ -1,137 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Robert Wagner 10/2002 -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHSigmabarTheta // -// // -// 2D-Histogram in Sigmabar and Theta // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHSigmabarTheta.h" - -#include - -#include - -#include "MMcEvt.hxx" -#include "MSigmabar.h" - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHSigmabarTheta); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets name and title of the histogram. -// -MHSigmabarTheta::MHSigmabarTheta(const char *name, const char *title) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHSigmabarTheta"; - fTitle = title ? title : "3-D histogram in sigmabar and theta"; - - fHist.SetDirectory(NULL); - - fHist.SetTitle("3D-plot of sigmabar and theta"); - fHist.SetXTitle("\\theta [\\circ]"); - fHist.SetYTitle("\\bar{\\sigma}"); - fHist.SetZTitle("N"); -} - -// -------------------------------------------------------------------------- -// -// Set binnings and prepare filling of the histogram -// -Bool_t MHSigmabarTheta::SetupFill(const MParList *plist) -{ - fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << "MMcEvt not found... aborting." << endl; - return kFALSE; - } - - fSigmabar = (MSigmabar*)plist->FindObject("MSigmabar"); - if (!fSigmabar) - { - *fLog << err << "MSigmabar not found... aborting." << endl; - return kFALSE; - } - - MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); - if (!binstheta) - { - *fLog << err << "BinningTheta [MBinning] not found... aborting." << endl; - return kFALSE; - } - MBinning* binssigmabar = (MBinning*)plist->FindObject("BinningSigmabar"); - if (!binssigmabar) - { - *fLog << err << "BinningSigmabar [MBinning] not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHist, binstheta, binssigmabar); - - fHist.Sumw2(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the histogram -// -Bool_t MHSigmabarTheta::Fill(const MParContainer *par, const Stat_t w) -{ - fHist.Fill(fMcEvt->GetTheta()*kRad2Deg, fSigmabar->GetSigmabar(), w); - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHSigmabarTheta::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - fHist.Draw(opt); - - pad->Modified(); - pad->Update(); -} - Index: trunk/MagicSoft/Mars/mhist/MHSigmabarTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmabarTheta.h (revision 6973) +++ (revision ) @@ -1,47 +1,0 @@ -#ifndef MARS_MHSigmabarTheta -#define MARS_MHSigmabarTheta - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH2 -#include -#endif - -class MMcEvt; -class MSigmabar; -class TH2D; -class MParList; - - -class MHSigmabarTheta : public MH -{ -private: - MMcEvt *fMcEvt; //! - MSigmabar *fSigmabar; //! - - TH2D fHist; - -public: - MHSigmabarTheta(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - const TH2D *GetHist() { return &fHist; } - const TH2D *GetHist() const { return &fHist; } - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void Draw(Option_t *option=""); - - ClassDef(MHSigmabarTheta, 1) //3D-histogram in alpha, Energy and theta -}; - -#endif - - - - - Index: trunk/MagicSoft/Mars/mhist/MHSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSupercuts.cc (revision 6973) +++ (revision ) @@ -1,225 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -/////////////////////////////////////////////////////////////////////// -// -// MHSupercuts -// -// This class contains histograms for the supercuts -// -// the histograms are filled during the optimization of the supercuts -// -/////////////////////////////////////////////////////////////////////// -#include "MHSupercuts.h" - -#include - -#include -#include -#include -#include - -#include "MParList.h" -#include "MHFindSignificance.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHSupercuts); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Setup four histograms for Alpha, and Dist -// -MHSupercuts::MHSupercuts(const char *name, const char *title) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHSupercuts"; - fTitle = title ? title : "Container for histograms for the supercuts"; - - - fDegree = new TH1F("Degree", "Degree of polynomial", 5, -0.5, 4.5); - fProb = new TH1F("Prob", "chi2 probability", 40, 0, 1.0); - fNdf = new TH1F("NDF", "NDF of polynomial fit", 60, -0.5, 59.5); - fGamma = new TH1F("gamma", "gamma", 40, 0.0, 8.0); - fNexNon = new TH1F("NexNon", "Nex / Non", 50, 0.0, 1.0); - fSigLiMa= new TH1F("Significance", "significance of gamma signal", - 60, 0.0, 120.0); - fSigtoBackg= new TH2F("SigtoBackg", "Significance vs signal/backg ratio", - 50, 0.0, 10.0, 60, 0.0, 120.0); - fSigDegree = new TH2F("SigDegree", "Significance vs Degree of polynomial", - 5, -0.5, 4.5, 60, 0.0, 120.0); - fSigNbins = new TH2F("SigNbins", "Significance vs number of bins", - 40, -0.5, 79.5, 60, 0.0, 120.0); - - fDegree->SetDirectory(NULL); - fProb->SetDirectory(NULL); - fNdf->SetDirectory(NULL); - fGamma->SetDirectory(NULL); - fNexNon->SetDirectory(NULL); - fSigLiMa->SetDirectory(NULL); - fSigtoBackg->SetDirectory(NULL); - fSigDegree->SetDirectory(NULL); - fSigNbins->SetDirectory(NULL); - - fDegree->SetXTitle("order of polynomial"); - fProb->SetXTitle("chi2 probability of polynomial fit"); - fNdf->SetXTitle("NDF of polynomial fit"); - fGamma->SetXTitle("gamma"); - fNexNon->SetXTitle("Nex / Non"); - fSigLiMa->SetXTitle("significance"); - - fSigtoBackg->SetXTitle("signa./background ratio"); - fSigtoBackg->SetYTitle("significance"); - - fSigDegree->SetXTitle("order of polynomial"); - fSigDegree->SetYTitle("significance"); - - fSigNbins->SetXTitle("number of bins"); - fSigNbins->SetYTitle("significance"); - - fDegree->SetYTitle("Counts"); - fProb->SetYTitle("Counts"); - fNdf->SetYTitle("Counts"); - fGamma->SetYTitle("Counts"); - fNexNon->SetYTitle("Counts"); - fSigLiMa->SetYTitle("Counts"); - - fSigtoBackg->SetZTitle("Counts"); - fSigDegree->SetZTitle("Counts"); - fSigNbins->SetZTitle("Counts"); -} - -// -------------------------------------------------------------------------- -// -// Delete the histograms -// -MHSupercuts::~MHSupercuts() -{ - delete fDegree; - delete fProb; - delete fNdf; - delete fGamma; - delete fNexNon; - delete fSigLiMa; - - delete fSigtoBackg; - delete fSigDegree; - delete fSigNbins; -} - -// -------------------------------------------------------------------------- -// -// Fill the histograms from the 'MHFindSignificance' container -// -Bool_t MHSupercuts::Fill(const MParContainer *par, const Stat_t w) -{ - if (!par) - { - *fLog << err << "MHSupercuts::Fill: Pointer (!=NULL) expected." << endl; - return kFALSE; - } - - MHFindSignificance &h = *(MHFindSignificance*)par; - - fDegree ->Fill(h.GetDegree( ), w); - fProb ->Fill(h.GetProb(), w); - fNdf ->Fill(h.GetNdf(), w); - fGamma ->Fill(h.GetGamma(), w); - - Double_t ratio = h.GetNon()>0.0 ? h.GetNex()/h.GetNon() : 0.0; - fNexNon ->Fill(ratio, w); - - fSigLiMa ->Fill(h.GetSigLiMa(), w); - - Double_t sigtobackg = h.GetNbg()!=0.0 ? h.GetNex() / h.GetNbg() : 0.0; - fSigtoBackg->Fill(sigtobackg, h.GetSigLiMa(), w); - - fSigDegree ->Fill(h.GetDegree(), h.GetSigLiMa(), w); - fSigNbins ->Fill(h.GetMbins(), h.GetSigLiMa(), w); - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Creates a new canvas and draws the two histograms into it. -// Be careful: The histograms belongs to this object and won't get deleted -// together with the canvas. -// -void MHSupercuts::Draw(Option_t *) -{ - //TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - //pad->SetBorderMode(0); - //AppendPad(""); - - TCanvas *pad = new TCanvas("Supercuts", "Supercut plots", 900, 900); - gROOT->SetSelectedPad(NULL); - - pad->Divide(3, 3); - - pad->cd(1); - gPad->SetBorderMode(0); - fDegree->Draw(); - - pad->cd(2); - gPad->SetBorderMode(0); - fProb->Draw(); - - pad->cd(3); - gPad->SetBorderMode(0); - fNdf->Draw(); - - pad->cd(4); - gPad->SetBorderMode(0); - fGamma->Draw(); - - pad->cd(5); - gPad->SetBorderMode(0); - fNexNon->Draw(); - - pad->cd(6); - gPad->SetBorderMode(0); - fSigLiMa->Draw(); - - pad->cd(7); - gPad->SetBorderMode(0); - fSigtoBackg->Draw(); - - pad->cd(8); - gPad->SetBorderMode(0); - fSigDegree->Draw(); - - pad->cd(9); - gPad->SetBorderMode(0); - fSigNbins->Draw(); - - pad->Modified(); - pad->Update(); -} Index: trunk/MagicSoft/Mars/mhist/MHSupercuts.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSupercuts.h (revision 6973) +++ (revision ) @@ -1,42 +1,0 @@ -#ifndef MARS_MHSupercuts -#define MARS_MHSupercuts - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH1F; -class TH2F; - -class MHSupercuts : public MH -{ -private: - TH1F *fDegree; // order of polynomial for background fit - TH1F *fProb; // chi2 probability of polynomial fit - TH1F *fNdf; // NDF of polynomial fit - TH1F *fGamma; // Nbg = gamma * Noff - TH1F *fNexNon; // no.of excess events / no.of events in signal region - TH1F *fSigLiMa; // significance of gamma signal - - TH2F *fSigtoBackg; // significance vs signal to background ratio (Nex/Nbg) - TH2F *fSigDegree; // significance vs order of polynomial - TH2F *fSigNbins; // significance vs number of bins - - -public: - MHSupercuts(const char *name=NULL, const char *title=NULL); - ~MHSupercuts(); - - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - void Draw(Option_t *opt=NULL); - - ClassDef(MHSupercuts, 1) // Container which holds histograms for the supercuts -}; - -#endif - - - - - Index: trunk/MagicSoft/Mars/mhist/MHThetabarTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHThetabarTheta.cc (revision 6973) +++ (revision ) @@ -1,121 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 1/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHThetabarTheta // -// // -// calculates the average Theta for different bins in Theta // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHThetabarTheta.h" - -#include - -#include "MMcEvt.hxx" - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHThetabarTheta); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets name and title of the histogram. -// -MHThetabarTheta::MHThetabarTheta(const char *name, const char *title) - : fHist() -{ - // - // set the name and title of this object - // - fName = name ? name : "MHThetabarTheta"; - fTitle = title ? title : "1-D profile histogram Thetabar vs. Theta"; - - fHist.SetDirectory(NULL); - - fHist.SetXTitle("\\Theta [\\circ]"); - fHist.SetYTitle("\\bar{\\Theta} [ \\circ]"); -} - -// -------------------------------------------------------------------------- -// -// Set the binnings and prepare the filling of the histogram -// -Bool_t MHThetabarTheta::SetupFill(const MParList *plist) -{ - fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; - return kFALSE; - } - - const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); - if (!binstheta ) - { - *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHist, binstheta); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHThetabarTheta::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - fHist.Draw(opt); - - pad->Modified(); - pad->Update(); -} - -// -------------------------------------------------------------------------- -// -// Fill the histogram -// -Bool_t MHThetabarTheta::Fill(const MParContainer *par, const Stat_t w) -{ - const Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg; - - fHist.Fill(theta, theta, w); - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mhist/MHThetabarTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHThetabarTheta.h (revision 6973) +++ (revision ) @@ -1,52 +1,0 @@ -#ifndef MARS_MHThetabarTheta -#define MARS_MHThetabarTheta - -#ifndef MARS_MH -#include "MH.h" -#endif -#ifndef ROOT_TProfile -#include -#endif - -class MTime; -class MMcEvt; - -class MParList; - -class MHThetabarTheta : public MH -{ -private: - MTime *fTime; //! - MMcEvt *fMcEvt; //! - - TProfile fHist; - -public: - MHThetabarTheta(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - const TProfile *GetHist() { return &fHist; } - const TProfile *GetHist() const { return &fHist; } - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void Draw(Option_t *option=""); - - ClassDef(MHThetabarTheta, 0) //Profile histogram Thetabar vs. time - -}; - -#endif - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mhist/MHThetabarTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHThetabarTime.cc (revision 6973) +++ (revision ) @@ -1,127 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 3/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHThetabarTime // -// // -// calculates the average Theta for different bins in time // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHThetabarTime.h" - -#include - -#include "MTime.h" -#include "MMcEvt.hxx" - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHThetabarTime); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets name and title of the histogram. -// -MHThetabarTime::MHThetabarTime(const char *name, const char *title) - : fHist() -{ - // - // set the name and title of this object - // - fName = name ? name : "MHThetabarTime"; - fTitle = title ? title : "1-D profile histogram Thetabar vs. time"; - - fHist.SetDirectory(NULL); - - fHist.SetXTitle("t [s]"); - fHist.SetYTitle("\\bar{\\Theta} [ \\circ]"); -} - -// -------------------------------------------------------------------------- -// -// Set the binnings and prepare the filling of the histogram -// -Bool_t MHThetabarTime::SetupFill(const MParList *plist) -{ - fTime = (MTime*)plist->FindObject("MTime"); - if (!fTime) - { - *fLog << err << dbginf << "MTime not found... aborting." << endl; - return kFALSE; - } - - fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; - return kFALSE; - } - - const MBinning* binstime = (MBinning*)plist->FindObject("BinningTime"); - if (!binstime ) - { - *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHist, binstime); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHThetabarTime::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - fHist.Draw(opt); - - pad->Modified(); - pad->Update(); -} - -// -------------------------------------------------------------------------- -// -// Fill the histogram -// -Bool_t MHThetabarTime::Fill(const MParContainer *par, const Stat_t w) -{ - fHist.Fill(*fTime, fMcEvt->GetTheta()*kRad2Deg, w); - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mhist/MHThetabarTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHThetabarTime.h (revision 6973) +++ (revision ) @@ -1,44 +1,0 @@ -#ifndef MARS_MHThetabarTime -#define MARS_MHThetabarTime - -#ifndef MARS_MH -#include "MH.h" -#endif -#ifndef ROOT_TProfile -#include "TProfile.h" -#endif -#ifndef ROOT_TH2 -#include "TH2.h" -#endif - -class MTime; -class MMcEvt; - -class MParList; - -class MHThetabarTime : public MH -{ -private: - MTime *fTime; //! - MMcEvt *fMcEvt; - - TProfile fHist; - -public: - MHThetabarTime(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - const TProfile *GetHist() { return &fHist; } - const TProfile *GetHist() const { return &fHist; } - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void Draw(Option_t *option=""); - - ClassDef(MHThetabarTime, 0) //Profile histogram Thetabar vs. time - -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/MHTrigLvl0.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHTrigLvl0.cc (revision 6973) +++ (revision ) @@ -1,234 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Abelardo Moralejo 06/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MHTrigLvl0 -// -// This is intended to be a sort of "level 0 trigger display". What it really -// does is to store the number of events of a data file in which each pixel -// has gone above a given threshold (fPixelThreshold) which is chosen when -// calling the constructor. Displaying a camera view with these values can -// help identify noisy pixels. See the macro pixfixrate.C to see an example -// of its use. Many things are to be fixed. Only inner pixels are shown now -// (which are anyhow those involved in the trigger), and the camera geometry -// (number of pixels, and how many inner ones) is not yet read from the input -// file. -// The "pedestal" we are using is just the signal in the first ADC slice -// (seems reasonable from the inspection of the available test data files). -// -// -///////////////////////////////////////////////////////////////////////////// -#include "MHTrigLvl0.h" - -#include - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" -#include "MRawEvtData.h" -#include "MRawEvtPixelIter.h" -#include "MCamDisplay.h" - -#include "MGeomCam.h" -#include "MGeomPix.h" - -ClassImp(MHTrigLvl0); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Reset all pixels to 0 and reset fEntries to 0. -// -void MHTrigLvl0::Clear() -{ - fSum.Reset(); - for (int i=0; i<577; i++) - { - fSum.AddPixel(i, 0, 0); - fSum[i].SetPixelUnused(); - } - fSum.FixSize(); - fEntries = 0; -} - -// -------------------------------------------------------------------------- -// -// Initialize the name and title of the task. -// Resets the sum histogram -// -MHTrigLvl0::MHTrigLvl0(const Float_t pixelthreshold, - const char *name, const char *title) - : fCam(NULL), fRawEvt(NULL), fDispl(NULL) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHTrigLvl0"; - fTitle = title ? title : "Number of hits above per pixel"; - fPixelThreshold = pixelthreshold; - - Clear(); -} - -// -------------------------------------------------------------------------- -// -// Delete the corresponding camera display if available -// -MHTrigLvl0::~MHTrigLvl0() -{ - if (fDispl) - delete fDispl; -} - -// -------------------------------------------------------------------------- -// -// Get the event (MRawEvtData) the histogram might be filled with. If -// it is not given, it is assumed, that it is filled with the argument -// of the Fill function. -// Looks for the camera geometry MGeomCam and resets the sum histogram. -// -Bool_t MHTrigLvl0::SetupFill(const MParList *plist) -{ - fRawEvt = (MRawEvtData*)plist->FindObject("MRawEvtData"); - if (!fRawEvt) - { - *fLog << dbginf << "MRawEvtData not found... aborting." << endl; - return kFALSE; - } - - fCam = (MGeomCam*)plist->FindObject("MGeomCam"); - if (!fCam) - *fLog << warn << GetDescriptor() << ": No MGeomCam found." << endl; - - Clear(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the histograms with data from a MCerPhotEvt-Container. -// -Bool_t MHTrigLvl0::Fill(const MParContainer *par, const Stat_t w) -{ - MRawEvtData *rawevt = par ? (MRawEvtData*) par : fRawEvt; - if (!rawevt) - { - *fLog << err << dbginf << "No MRawEvtData found..." << endl; - return kFALSE; - } - - MRawEvtPixelIter pixel(rawevt); - - while(pixel.Next()) - { - const UInt_t pixid = pixel.GetPixelId(); - - // FIXME: number of inner pixels should be read from file - if (pixid > 396) - break; - if (pixid == 0) - continue; - - fSum[pixid].SetPixelUsed(); - - // - // FIXME: we now use as "pedestal" the value of the first ADC slice! - // - Float_t baseline = rawevt->GetNumHiGainSamples() * - pixel.GetHiGainSamples()[0]; - - Float_t pixel_signal = pixel.GetSumHiGainSamples() - baseline; - - Float_t pixel_is_on = ( pixel_signal > fPixelThreshold)? 1. : 0.; - - fSum[pixid].AddNumPhotons(pixel_is_on); - } - - fEntries++; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Set to Unused outer pixels. -// -Bool_t MHTrigLvl0::Finalize() -{ - // - // Show only pixels in the inner region of the camera: - // (otherwise, problem with the too different ranges) - // - for (Int_t i=0; i<577; i++) - // FIXME: read number of total and inner pixels from file - { - if (i > 396) - fSum[i].SetPixelUnused(); - } - - // fSum.Scale(fEntries); Now disabled, scale was not readable otherwise. - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Draw the present 'fill status' -// -void MHTrigLvl0::Draw(Option_t *) -{ - if (!fCam) - { - *fLog << warn << "WARNING - Cannot draw " << GetDescriptor() << ": No Camera Geometry available." << endl; - return; - } - - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this, 750, 600); - pad->SetBorderMode(0); - - AppendPad(""); -} - -// -------------------------------------------------------------------------- -// -// If a camera display is not yet assigned, assign a new one. -// -void MHTrigLvl0::Paint(Option_t *option) -{ - if (!fCam) - { - *fLog << warn << "WARNING - Cannot paint " << GetDescriptor() << ": No Camera Geometry available." << endl; - return; - } - - if (!fDispl) - fDispl = new MCamDisplay(fCam); - - fDispl->Fill(fSum); - fDispl->Paint(); -} Index: trunk/MagicSoft/Mars/mhist/MHTrigLvl0.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHTrigLvl0.h (revision 6973) +++ (revision ) @@ -1,46 +1,0 @@ -#ifndef MARS_MHTrigLvl0 -#define MARS_MHTrigLvl0 - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef MARS_MCerPhotEvt -#include "MCerPhotEvt.h" -#endif - -class TH1D; -class MCamDisplay; -class MRawEvtData; - -class MHTrigLvl0 : public MH -{ -private: - MCerPhotEvt fSum; // storing the sum of triggers - Int_t fEntries; // number of entries in the histogram - MGeomCam *fCam; // the present geometry - MRawEvtData *fRawEvt; //! ADC info of the current event - MCamDisplay *fDispl; //! the camera display - Float_t fPixelThreshold; // Threshold (ADC counts a.p.) to consider - // a pixel as "fired". - -public: - MHTrigLvl0(const Float_t pixelthreshold = 0, - const char *name=NULL, const char *title=NULL); - ~MHTrigLvl0(); - - void Clear(); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - Bool_t Finalize(); - - const MCerPhotEvt &GetSum() const { return fSum; } - - void Draw(Option_t *opt=""); - void Paint(Option_t *option=""); - - ClassDef(MHTrigLvl0, 1) // Histogram to sum level 0 triggers in all pixels -}; - -#endif