Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2295) +++ trunk/MagicSoft/Mars/Changelog (revision 2296) @@ -111,6 +111,6 @@ - add member function GetSignificance() - * mhist/MHMatrix.cc - - add MProgressBar in Fill() + * mhist/MHMatrix.cc + - add MProgressBar in Fill() * mmontecarlo/MMcEnergyEst.h Index: trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc (revision 2295) +++ trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc (revision 2296) @@ -172,4 +172,5 @@ for (UShort_t i=0; iGetErrorPhot(); // FIXME: perr[i] ??? + num++; } - num++; } Index: trunk/MagicSoft/Mars/mbase/MArgs.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MArgs.cc (revision 2295) +++ trunk/MagicSoft/Mars/mbase/MArgs.cc (revision 2296) @@ -38,4 +38,6 @@ ClassImp(MArgsEntry); ClassImp(MArgs); + +using namespace std; void MArgsEntry::Print(const Option_t *o) const Index: trunk/MagicSoft/Mars/mbase/MDirIter.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MDirIter.h (revision 2295) +++ trunk/MagicSoft/Mars/mbase/MDirIter.h (revision 2296) @@ -38,8 +38,8 @@ AddDirectory(o->GetName(), o->GetTitle()); } - MDirIter(const char *dir, const char *filter="") : fNext(&fList), fDirPtr(NULL) + MDirIter(const char *dir, const char *filter="", Int_t rec=0) : fNext(&fList), fDirPtr(NULL) { fList.SetOwner(); - AddDirectory(dir, filter); + AddDirectory(dir, filter, rec); } ~MDirIter() Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 2295) +++ trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 2296) @@ -125,5 +125,5 @@ // Add this file as the last entry in the chain // -void MCT1ReadPreProc::AddFile(const char *txt) +Int_t MCT1ReadPreProc::AddFile(const char *txt, int) { const char *name = gSystem->ExpandPathName(txt); @@ -135,5 +135,5 @@ { *fLog << warn << "WARNING - Problem reading header... ignored." << endl; - return; + return 0; } @@ -142,5 +142,5 @@ { *fLog << warn << "WARNING - File contains no data... ignored." << endl; - return; + return 0; } @@ -150,4 +150,5 @@ fFileNames->AddLast(new TNamed(txt, "")); + return 1; } Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 2295) +++ trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 2296) @@ -92,5 +92,5 @@ ~MCT1ReadPreProc(); - void AddFile(const char *fname); + Int_t AddFile(const char *fname, int i=0); UInt_t GetEntries() { return fEntries; } Index: trunk/MagicSoft/Mars/mfileio/MRead.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MRead.cc (revision 2295) +++ trunk/MagicSoft/Mars/mfileio/MRead.cc (revision 2296) @@ -39,4 +39,6 @@ #include "MLog.h" #include "MLogManip.h" + +#include "MDirIter.h" ClassImp(MRead); @@ -98,2 +100,16 @@ return i!=0; } + +Int_t MRead::AddFiles(MDirIter &next) +{ + Int_t rc=0; + TString n; + while (!(n=next()).IsNull()) + { + const Int_t num = AddFile(n); + if (num>0) + rc += num; + } + return rc; +} + Index: trunk/MagicSoft/Mars/mfileio/MRead.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/MRead.h (revision 2295) +++ trunk/MagicSoft/Mars/mfileio/MRead.h (revision 2296) @@ -7,4 +7,5 @@ class MFilter; +class MDirIter; class MRead : public MTask @@ -22,5 +23,6 @@ MFilter *GetSelector() { return fSelector; } - Int_t AddFile(const char *fname, Int_t entries=-1) { return -1; } + virtual Int_t AddFile(const char *fname, Int_t entries=-1) { return -1; } + Int_t AddFiles(MDirIter &next); Bool_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); Index: trunk/MagicSoft/Mars/mfileio/MReadRflFile.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MReadRflFile.cc (revision 2295) +++ trunk/MagicSoft/Mars/mfileio/MReadRflFile.cc (revision 2296) @@ -399,5 +399,5 @@ // Add this file as the last entry in the chain // -void MReadRflFile::AddFile(const char *txt) +Int_t MReadRflFile::AddFile(const char *txt, int) { const char *name = gSystem->ExpandPathName(txt); @@ -424,4 +424,5 @@ */ fFileNames->AddLast(new TNamed(txt, "")); + return 1; } Index: trunk/MagicSoft/Mars/mfileio/MReadRflFile.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/MReadRflFile.h (revision 2295) +++ trunk/MagicSoft/Mars/mfileio/MReadRflFile.h (revision 2296) @@ -48,5 +48,5 @@ ~MReadRflFile(); - void AddFile(const char *fname); + Int_t AddFile(const char *fname, int i=0); Bool_t Rewind() { fNumFile=0; return kTRUE; } Index: trunk/MagicSoft/Mars/mfileio/MReadTree.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 2295) +++ trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 2296) @@ -121,6 +121,7 @@ if (fname) - if (fChain->Add(fname)>0) - SetBit(kChainWasChanged); + AddFile(fname); +// if (fChain->Add(fname)>0) +// SetBit(kChainWasChanged); } @@ -586,5 +587,5 @@ if (!GetEntries()) { - *fLog << warn << dbginf << "No entries found in file(s)" << endl; + *fLog << warn << GetDescriptor() << ": No entries found in file(s)" << endl; return kFALSE; } Index: trunk/MagicSoft/Mars/mhist/MH.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MH.cc (revision 2295) +++ trunk/MagicSoft/Mars/mhist/MH.cc (revision 2296) @@ -796,10 +796,12 @@ // Otherwise the present gPad is returned. // -TVirtualPad *MH::GetNewPad(Option_t *opt) -{ - TString str(opt); - - if (!str.Contains("nonew", TString::kIgnoreCase)) +TVirtualPad *MH::GetNewPad(TString &opt) +{ + opt.ToLower(); + + if (!opt.Contains("nonew")) return NULL; + + opt.ReplaceAll("nonew", ""); return gPad; @@ -813,9 +815,8 @@ TObject *MH::Clone(const char *name) const { - Bool_t store = TH1::AddDirectoryStatus(); + const Bool_t store = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); - TObject *o = MParContainer::Clone(name); - TH1::AddDirectory(store); @@ -831,5 +832,7 @@ TObject *MH::DrawClone(Option_t *opt, Int_t w, Int_t h) const { - TVirtualPad *p = GetNewPad(opt); + TString option(opt); + + TVirtualPad *p = GetNewPad(option); if (!p) p = MakeDefCanvas(this, w, h); @@ -839,5 +842,5 @@ gROOT->SetSelectedPad(NULL); - TObject *o = MParContainer::DrawClone(opt); + TObject *o = MParContainer::DrawClone(option); o->SetBit(kCanDelete); return o; Index: trunk/MagicSoft/Mars/mhist/MH.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MH.h (revision 2295) +++ trunk/MagicSoft/Mars/mhist/MH.h (revision 2296) @@ -75,5 +75,5 @@ } - static TVirtualPad *GetNewPad(Option_t *opt); + static TVirtualPad *GetNewPad(TString &opt); static void FindGoodLimits(Int_t nbins, Int_t &newbins, Double_t &xmin, Double_t &xmax, Bool_t isInteger); Index: trunk/MagicSoft/Mars/mhist/MHHadronness.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHHadronness.cc (revision 2295) +++ trunk/MagicSoft/Mars/mhist/MHHadronness.cc (revision 2296) @@ -430,4 +430,17 @@ // -------------------------------------------------------------------------- // +// Search the hadronness corresponding to a given hadron acceptance. +// +Double_t MHHadronness::GetHadronness(Double_t acchad) const +{ + for (int i=1; iGetNbinsX()+1; i++) + if (fIntPhness->GetBinContent(i)>acchad) + return fIntPhness->GetBinLowEdge(i); + + return -1; +} + +// -------------------------------------------------------------------------- +// // Print the corresponding Gammas Acceptance for a hadron acceptance of // 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum @@ -448,12 +461,12 @@ *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl; - *fLog << "acc(hadron) acc(gamma) acc(g)/acc(h)" << endl <Remove(&tlist); @@ -1097,19 +1104,33 @@ void MHMatrix::ShuffleRows(UInt_t seed) { - TRandom rnd(seed); - - for (Int_t irow = 0; irow < fNumRow; irow++) - { - Int_t jrow = rnd.Integer(fNumRow); - - for (Int_t icol = 0; icol < fM.GetNcols(); icol++) - { - Real_t tmp = fM(irow,icol); - fM(irow,icol) = fM(jrow,icol); - fM(jrow,icol) = tmp; - } - } - - *fLog << warn << this->GetName() << " : Attention! Matrix rows have been shuffled." << endl; - -} + TRandom rnd(seed); + + TVector v(fM.GetNcols()); + TVector tmp(fM.GetNcols()); + for (Int_t irow = 0; irow=fM.GetNrows()) + { + *fLog << warn << GetDescriptor() << ": Warning - " << num << " >= rows=" << fM.GetNrows() << endl; + return; + } + + const TMatrix m(fM); + fM.ResizeTo(num, fM.GetNcols()); + + TVector tmp(fM.GetNcols()); + for (UInt_t irow=0; irow +! Author(s): Harald Kornmayer, 1/2001 +! Author(s): Nadia Tonello, 4/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MImgCleanTGB +// +// The Image Cleaning task selects the pixels you use for the Hillas +// parameters calculation. +// +// There are two methods to make the selection: the standard one, as done +// in the analysis of CT1 data, and the democratic one, as suggested by +// W.Wittek. The number of photo-electrons of a pixel is compared with the +// pedestal RMS of the pixel itself (standard method) or with the average +// RMS of the inner pixels (democratic method). +// In both cases, the possibility to have a camera with pixels of +// different area is taken into account. +// The too noisy pixels can be recognized and eventally switched off +// (Unmap: set blind pixels to UNUSED) separately, using the +// MBlindPixelCalc Class. In the MBlindPixelCalc class there is also the +// function to replace the value of the noisy pixels with the interpolation +// of the content of the neighbors (SetUseInterpolation). +// +// Example: +// ... +// MBlindPixelCalc blind; +// blind.SetUseInterpolation(); +// blind.SetUseBlindPixels(); +// +// MImgCleanTGB clean; +// ... +// tlist.AddToList(&blind); +// tlist.AddToList(&clean); +// +// Look at the MBlindPixelCalc Class for more details. +// +// Starting point: default values ---------------------------------------- +// +// When an event is read, before the image cleaning, all the pixels that +// are in MCerPhotEvt are set as USED and NOT CORE. All the pixels belong +// to RING number 1 (like USED pixels). +// Look at MCerPhotPix.h to see how these informations of the pixel are +// stored. +// The default cleaning METHOD is the STANDARD one and the number of the +// rings around the CORE pixel it analyzes is 1. Look at the Constructor +// of the class in MImgCleanTGB.cc to see (or change) the default values. +// +// Example: To modify this setting, use the member functions +// SetMethod(MImgCleanTGB::kDemocratic) and SetCleanRings(UShort_t n). +// +// MImgCleanTGB:CleanStep1 ----------------------------------------------- +// +// The first step of cleaning defines the CORE pixels. The CORE pixels are +// the ones which contain the informations about the core of the electro- +// magnetic shower. +// The ratio (A_0/A_i) is calculated from fCam->GetPixRatio(i). A_0 is +// the area of the central pixel of the camera, A_i is the area of the +// examined pixel. In this way, if we have a MAGIC-like camera, with the +// outer pixels bigger than the inner ones, the level of cleaning in the +// two different regions is weighted. +// This avoids problems of deformations of the shower images. +// The signal S_i and the pedestal RMS Prms_i of the pixel are called from +// the object MCerPhotPix. +// If (default method = kStandard) +//Begin_Html +//       +//End_Html +// the pixel is set as CORE pixel. L_1 (n=1) is called "first level of +// cleaning" (default: fCleanLvl1 = 3). +// All the other pixels are set as UNUSED and belong to RING 0. +// After this point, only the CORE pixels are set as USED, with RING +// number 1. +// +// MImgCleanTGB:CleanStep2 ---------------------------------------------- +// +// The second step of cleaning looks at the isolated CORE pixels and sets +// them to UNUSED. An isolated pixel is a pixel without CORE neighbors. +// At the end of this point, we have set as USED only CORE pixels with at +// least one CORE neighbor. +// +// MImgCleanTGB:CleanStep3 ---------------------------------------------- +// +// The third step of cleaning looks at all the pixels (USED or UNUSED) that +// surround the USED pixels. +// If the content of the analyzed pixel survives at the second level of +// cleaning, i.e. if +//Begin_Html +//       +//End_Html +// the pixel is set as USED. L_2 (n=2) is called "second level of cleaning" +// (default:fCleanLvl2 = 2.5). +// +// When the number of RINGS to analyze is 1 (default value), only the +// pixels that have a neighbor CORE pixel are analyzed. +// +// There is the option to decide the number of times you want to repeat +// this procedure (number of RINGS analyzed around the core pixels = n). +// Every time the level of cleaning is the same (fCleanLvl2) and the pixel +// will belong to ring r+1, 1 < r < n+1. This is described in +// MImgCleanTGB:CleanStep4 . +// +// Dictionary and member functions --------------------------------------- +// +// Here there is the detailed description of the member functions and of +// the terms commonly used in the class. +// +// STANDARD CLEANING: +// ================= +// This is the method used for the CT1 data analysis. It is the default +// method of the class. +// The number of photo-electrons of a pixel (S_i) is compared to the +// pedestal RMS of the pixel itself (Prms_i). To have the comparison to +// the same photon density for all the pixels, taking into account they +// can have different areas, we have to keep in mind that the number of +// photons that hit each pixel, goes linearly with the area of the pixel. +// The fluctuations of the LONS are proportional to sqrt(A_i), so when we +// compare S_i with Prms_i, only a factor sqrt(A_0/A_i) is missing to +// have the same (N.photons/Area) threshold for all the pixels. +// +// !!WARNING: if noise independent from the +// pixel size (example: electronic noise) is introduced, +// then the noise fluctuations are no longer proportional +// to sqrt(A_i), and then the cut value (for a camera with +// pixels of different sizes) resulting from the above +// procedure would not be proportional to pixel size as we +// intend. In that case, democratic cleaning is preferred. +// +// If +//Begin_Html +//       +//End_Html +// the pixel survives the cleaning and it is set as CORE (when L_n is the +// first level of cleaning, fCleanLvl1) or USED (when L_n is the second +// level of cleaning, fCleanLvl2). +// +// Example: +// +// MImgCleanTGB clean; +// //creates a default Cleaning object, with default setting +// ... +// tlist.AddToList(&clean); +// // add the image cleaning to the main task list +// +// DEMOCRATIC CLEANING: +// =================== +// You use this cleaning method when you want to compare the number of +// photo-electons of each pixel with the average pedestal RMS +// (fInnerNoise = fSgb->GetSigmabarInner()) of the inner pixels (for the +// MAGIC camera they are the smaller ones): +//Begin_Html +//       +//End_Html +// In this case, the simple ratio (A_0/A_i) is used to weight the level of +// cleaning, because both the inner and the outer pixels (that in MAGIC +// have a different area) are compared to the same pedestal RMS, coming +// from the inner pixels. +// To calculate the average pedestal RMS of the inner pixels, you have to +// add to the main task list an object of type MSigmabarCalc before the +// MImgCleanTGB object. To know how the calculation of fInnerNoise is done +// look at the MSigmabarCalc Class. +// +// Example: +// +// MSigmabarCalc sbcalc; +// //creates an object that calcutates the average pedestal RMS +// MImgCleanTGB clean; +// ... +// tlist.AddToList(&sbcalc); +// tlist.AddToList(&clean); +// +// Member Function: SetMethod() +// ============================ +// When you call the MImgCleanTGB task, the default method is kStandard. +// +// If you want to switch to the kDemocratic method you have to +// call this member function. +// +// Example: +// +// MImgCleanTGB clean; +// //creates a default Cleaning object, with default setting +// +// clean.SetMethod(MImgCleanTGB::kDemocratic); +// //now the method of cleaning is changed to Democratic +// +// FIRST AND SECOND CLEANING LEVEL +// =============================== +// When you call the MImgCleanTGB task, the default cleaning levels are +// fCleanLvl1 = 3, fCleanLvl2 = 2.5. You can change them easily when you +// create the MImgCleanTGB object. +// +// Example: +// +// MImgCleanTGB clean(Float_t lvl1,Float_t lvl2); +// //creates a default cleaning object, but the cleaning levels are now +// //lvl1 and lvl2. +// +// RING NUMBER +// =========== +// The standard cleaning procedure is such that it looks for the +// informations of the boundary part of the shower only on the first +// neighbors of the CORE pixels. +// There is the possibility now to look not only at the firs neighbors +// (first ring),but also further away, around the CORE pixels. All the new +// pixels you can find with this method, are tested with the second level +// of cleaning and have to have at least an USED neighbor. +// +// They will be also set as USED and will be taken into account during the +// calculation of the image parameters. +// The only way to distinguish them from the other USED pixels, is the +// Ring number, that is bigger than 1. +// +// Example: You can decide how many rings you want to analyze using: +// +// MImgCleanTGB clean; +// //creates a default cleaning object (default number of rings =1) +// clean.SetCleanRings(UShort_t r); +// //now it looks r times around the CORE pixels to find new pixels with +// //signal. +// +// +// Input Containers: +// MGeomCam, MCerPhotEvt, MSigmabar +// +// Output Containers: +// MCerPhotEvt +// +///////////////////////////////////////////////////////////////////////////// +#include "MImgCleanTGB.h" + +#include // atof +#include // ofstream, SavePrimitive + +#include // TGFrame +#include // TGLabel +#include // TGTextEntry + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MSigmabar.h" + +#include "MGeomPix.h" +#include "MGeomCam.h" + +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" + +#include "MPedestalPix.h" +#include "MPedestalCam.h" + +#include "MGGroupFrame.h" // MGGroupFrame + +ClassImp(MImgCleanTGB); + +using namespace std; + +enum { + kImgCleanLvl1, + kImgCleanLvl2 +}; + +static const TString gsDefName = "MImgCleanTGB"; +static const TString gsDefTitle = "Task to perform image cleaning"; + +// -------------------------------------------------------------------------- +// +// Default constructor. Here you can specify the cleaning method and levels. +// If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma +// above mean) are used. +// Here you can also specify how many rings around the core pixels you want +// to analyze (with the fixed lvl2). The default value for "rings" is 1. +// +MImgCleanTGB::MImgCleanTGB(const Float_t lvl1, const Float_t lvl2, + const char *name, const char *title) + : fSgb(NULL), fCleaningMethod(kStandard), fCleanLvl1(lvl1), + fCleanLvl2(lvl2), fCleanRings(1) + +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + Print(); +} + + +Int_t MImgCleanTGB::CleanStep3b(MCerPhotPix &pix) +{ + const Int_t id = pix.GetPixId(); + + // + // check if the pixel's next neighbor is a core pixel. + // if it is a core pixel set pixel state to: used. + // + MGeomPix &gpix = (*fCam)[id]; + const Int_t nnmax = gpix.GetNumNeighbors(); + + Int_t rc = 0; + + for (Int_t j=0; jGetPixById(id2) && fEvt->IsPixelUsed(id2)) + rc++; + } + return rc; +} + +// -------------------------------------------------------------------------- +// +// Look for the boundary pixels around the core pixels +// if a pixel has more than 2.5 (clean level 2.5) sigma, and +// a core neigbor, it is declared as used. +// +void MImgCleanTGB::CleanStep3(Int_t num1, Int_t num2) +{ + const Int_t entries = fEvt->GetNumPixels(); + + Int_t *u = new Int_t[entries]; + + for (Int_t i=0; inum2) + pix.SetPixelUsed(); + } + + delete u; +} + +void MImgCleanTGB::CleanStep3(Byte_t *nb, Int_t num1, Int_t num2) +{ + const Int_t entries = fEvt->GetNumPixels(); + + for (Int_t i=0; inum2 && !pix.IsPixelUsed()) + { + const MGeomPix &gpix = (*fCam)[idx]; + const Int_t nnmax = gpix.GetNumNeighbors(); + for (Int_t j=0; jFindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + if (fCleaningMethod == kDemocratic) + { + fSgb = (MSigmabar*)pList->FindObject("MSigmabar"); + if (!fSgb) + { + *fLog << dbginf << "MSigmabar not found... aborting." << endl; + return kFALSE; + } + } + else + { + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << dbginf << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Cleans the image. +// +Int_t MImgCleanTGB::Process() +{ + const Int_t entries = fEvt->GetNumPixels(); + + Double_t sum = 0; + Double_t sq = 0; + Double_t w = 0; + Double_t w2 = 0; + for (Int_t i=0; iGetPixRatio(pix.GetPixId()); + const Float_t noise = (*fPed)[pix.GetPixId()].GetPedestalRms(); + + if (entry * TMath::Sqrt(factor) <= fCleanLvl2 * noise) + { + sum += entry*factor; + sq += entry*entry*factor*factor; + w += factor; + w2 += factor*factor; + } + } + + Double_t mean = sum/w; + Double_t sdev = sqrt(sq/w2 - mean*mean); + + Byte_t *nb = new Byte_t[1000]; + memset(nb, 0, 577); + + for (Int_t i=0; iGetPixRatio(idx); + + if (entry*factor > fCleanLvl1*sdev) + { + pix.SetPixelUsed(); + + const MGeomPix &gpix = (*fCam)[idx]; + const Int_t nnmax = gpix.GetNumNeighbors(); + for (Int_t j=0; j use with root >=3.02 <-- + * + + TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative); + TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative); + + */ + TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1); + TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2); + + // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight); + // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight); + + entry1->SetText("3.0"); + entry2->SetText("2.5"); + + entry1->Associate(f); + entry2->Associate(f); + + TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1"); + TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2"); + + l1->SetTextJustify(kTextLeft); + l2->SetTextJustify(kTextLeft); + + // + // Align the text of the label centered, left in the row + // with a left padding of 10 + // + TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10); + TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 5, 0, 10); + + // + // Add one entry field and the corresponding label to each line + // + f1->AddFrame(entry1); + f2->AddFrame(entry2); + + f1->AddFrame(l1, laylabel); + f2->AddFrame(l2, laylabel); + + f->AddFrame(f1, layframe); + f->AddFrame(f2, layframe); + + f->AddToList(entry1); + f->AddToList(entry2); + f->AddToList(l1); + f->AddToList(l2); + f->AddToList(laylabel); + f->AddToList(layframe); +} + +// -------------------------------------------------------------------------- +// +// Process the GUI Events comming from the two text entry fields. +// +Bool_t MImgCleanTGB::ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2) +{ + if (msg!=kC_TEXTENTRY || submsg!=kTE_ENTER) + return kTRUE; + + TGTextEntry *txt = (TGTextEntry*)FindWidget(param1); + + if (!txt) + return kTRUE; + + Float_t lvl = atof(txt->GetText()); + + switch (param1) + { + case kImgCleanLvl1: + fCleanLvl1 = lvl; + *fLog << "Cleaning level 1 set to " << lvl << " sigma." << endl; + return kTRUE; + + case kImgCleanLvl2: + fCleanLvl2 = lvl; + *fLog << "Cleaning level 2 set to " << lvl << " sigma." << endl; + return kTRUE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Implementation of SavePrimitive. Used to write the call to a constructor +// to a macro. In the original root implementation it is used to write +// gui elements to a macro-file. +// +void MImgCleanTGB::StreamPrimitive(ofstream &out) const +{ + out << " MImgCleanTGB " << GetUniqueName() << "("; + out << fCleanLvl1 << ", " << fCleanLvl2; + + if (fName!=gsDefName || fTitle!=gsDefTitle) + { + out << ", \"" << fName << "\""; + if (fTitle!=gsDefTitle) + out << ", \"" << fTitle << "\""; + } + out << ");" << endl; + + if (fCleaningMethod!=kDemocratic) + return; + + out << " " << GetUniqueName() << ".SetMethod(MImgCleanTGB::kDemocratic);" << endl; + + if (fCleanRings==1) + return; + + out << " " << GetUniqueName() << ".SetCleanRings(" << fCleanRings << ");" << endl; +} Index: trunk/MagicSoft/Mars/mimage/MImgCleanTGB.h =================================================================== --- trunk/MagicSoft/Mars/mimage/MImgCleanTGB.h (revision 2296) +++ trunk/MagicSoft/Mars/mimage/MImgCleanTGB.h (revision 2296) @@ -0,0 +1,66 @@ +#ifndef MARS_MImgCleanTGB +#define MARS_MImgCleanTGB + +#ifndef MARS_MGTask +#include "MGTask.h" +#endif + +class MGeomCam; +class MSigmabar; +class MCerPhotPix; +class MCerPhotEvt; +class MPedestalCam; + +class MGGroupFrame; + +class MImgCleanTGB : public MGTask +{ +public: + typedef enum { + kStandard, + kDemocratic + } CleaningMethod_t; + +private: + const MGeomCam *fCam; //! + MCerPhotEvt *fEvt; //! + MSigmabar *fSgb; //! + MPedestalCam *fPed; //! + + CleaningMethod_t fCleaningMethod; + + Float_t fCleanLvl1; + Float_t fCleanLvl2; + + UShort_t fCleanRings; + + Float_t fInnerNoise; //! + + void CreateGuiElements(MGGroupFrame *f); + void StreamPrimitive(ofstream &out) const; + + Int_t CleanStep3b(MCerPhotPix &pix); + void CleanStep3(Int_t num1, Int_t num2); + void CleanStep3(Byte_t *nb, Int_t num1, Int_t num2); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + +public: + MImgCleanTGB(const Float_t lvl1=3.0, const Float_t lvl2=2.5, + const char *name=NULL, const char *title=NULL); + void Print(Option_t *o="") const; + + Float_t GetCleanLvl1() const { return fCleanLvl1; } + Float_t GetCleanLvl2() const { return fCleanLvl2; } + UShort_t GetCleanRings() const { return fCleanRings;} + + void SetCleanRings(UShort_t r) { if(r==0) r=1; fCleanRings=r; } + void SetMethod(CleaningMethod_t m) { fCleaningMethod = m; } + + Bool_t ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2); + + ClassDef(MImgCleanTGB, 2) // task doing the image cleaning +}; + +#endif Index: trunk/MagicSoft/Mars/mimage/Makefile =================================================================== --- trunk/MagicSoft/Mars/mimage/Makefile (revision 2295) +++ trunk/MagicSoft/Mars/mimage/Makefile (revision 2296) @@ -29,4 +29,5 @@ SRCFILES = MImgCleanStd.cc \ + MImgCleanTGB.cc \ MCameraSmooth.cc \ MHillas.cc \ Index: trunk/MagicSoft/Mars/mmain/MStatusDisplay.cc =================================================================== --- trunk/MagicSoft/Mars/mmain/MStatusDisplay.cc (revision 2295) +++ trunk/MagicSoft/Mars/mmain/MStatusDisplay.cc (revision 2296) @@ -840,4 +840,7 @@ Bool_t MStatusDisplay::HasCanvas(const TCanvas *c) const { + if (!c) + return kFALSE; + if (gROOT->IsBatch()) return (Bool_t)fBatch->FindObject(c); Index: trunk/MagicSoft/Mars/mranforest/MHRanForest.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MHRanForest.cc (revision 2295) +++ trunk/MagicSoft/Mars/mranforest/MHRanForest.cc (revision 2296) @@ -68,5 +68,4 @@ fGraphSigma = new TGraph; fGraphSigma->SetTitle("Evolution of Standard deviation of estimated hadronness in tree combination"); - fGraphSigma->SetMaximum(1); fGraphSigma->SetMarkerStyle(kFullDotSmall); } @@ -113,4 +112,5 @@ { fNumEvent++; + Double_t hest=0; Double_t htrue=fMcEvt->GetPartId()==kGAMMA ? 0. : 1.; @@ -124,5 +124,7 @@ hest/=i+1; - fSigma[i]+=(htrue-hest)*(htrue-hest); + + const Double_t val = htrue-hest; + fSigma[i] += val*val; } @@ -145,11 +147,15 @@ for (Int_t i=0; iSetPoint(i,ip,ig); + fSigma[i] = TMath::Sqrt(fSigma[i]/fNumEvent); + + const Stat_t ig = fSigma[i]; + if (ig>max) max = ig; + if (igSetPoint(i, i+1, ig); } + + // This is used in root>3.04/? so that SetMaximum/Minimum can succeed + fGraphSigma->GetHistogram(); + fGraphSigma->SetMaximum(1.05*max); fGraphSigma->SetMinimum(0.95*min); @@ -180,5 +186,5 @@ return; - h->GetXaxis()->SetRangeUser(0, 1); + //h->GetXaxis()->SetRangeUser(0, fNumEvent+1); h->SetXTitle("No.of Trees"); h->SetYTitle("\\sigma of est.hadronness"); Index: trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc (revision 2295) +++ trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc (revision 2296) @@ -66,5 +66,4 @@ fGraphGini = new TGraph; fGraphGini->SetTitle("Importance of RF-input parameters measured by mean Gini decrease"); - fGraphGini->SetMaximum(1); fGraphGini->SetMarkerStyle(kFullDotSmall); } @@ -116,23 +115,27 @@ Bool_t MHRanForestGini::Finalize() { - Int_t n = fGini.GetSize(); + const Int_t n = fGini.GetSize(); fGraphGini->Set(n); - Stat_t max=0.; - Stat_t min=0.; + Stat_t max=0; + Stat_t min=0; for (Int_t i=0; iGetNumTrees(); - fGini[i]/=fRanForest->GetNumData(); + fGini[i] /= fRanForest->GetNumTrees(); + fGini[i] /= fRanForest->GetNumData(); - Stat_t ip = i+1.; - Stat_t ig = fGini[i]; + const Stat_t ip = i+1; + const Stat_t ig = fGini[i]; - max=TMath::Max(max,ig); - min=TMath::Min(min,ig); + if (ig>max) max=ig; + if (igSetPoint(i,ip,ig); } + + // This is used in root>3.04/? so that SetMaximum/Minimum can succeed + fGraphGini->GetHistogram(); + fGraphGini->SetMaximum(1.05*max); fGraphGini->SetMinimum(0.95*min); @@ -163,5 +166,5 @@ return; - h->GetXaxis()->SetRangeUser(0, 1); + //h->GetXaxis()->SetRangeUser(0, fGini.GetSize()+1); h->SetXTitle("No.of RF-input parameter"); h->SetYTitle("Mean decrease in Gini-index [a.u.]"); Index: trunk/MagicSoft/Mars/mranforest/MRanForest.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForest.cc (revision 2295) +++ trunk/MagicSoft/Mars/mranforest/MRanForest.cc (revision 2296) @@ -103,8 +103,8 @@ Int_t ntree=0; - TIter forest(fForest); + TIter Next(fForest); MRanTree *tree; - while ((tree=(MRanTree*)forest.Next())) + while ((tree=(MRanTree*)Next())) { fTreeHad[ntree]=tree->TreeHad(event); @@ -112,5 +112,5 @@ ntree++; } - return hadroness/Double_t(ntree); + return hadroness/ntree; } @@ -190,29 +190,23 @@ Bool_t MRanForest::GrowForest() { - Int_t ninbag=0; - TArrayI datsortinbag; - TArrayF classpopw; - TArrayI jinbag; - TArrayF winbag; - - jinbag.Set(fNumData); - winbag.Set(fNumData); - classpopw.Set(2); - - TMatrix hadrons=fHadrons->GetM(); - TMatrix gammas=fGammas->GetM(); + if(!gRandom) + { + *fLog << err << dbginf << "gRandom not initialized... aborting." << endl; + return kFALSE; + } fTreeNo++; // initialize running output - if(fTreeNo==1) - { - cout< overestim. of error)"< overestim. of error)" << endl; + } + + TArrayF classpopw(2); + TArrayI jinbag(fNumData); // Initialization includes filling with 0 + TArrayF winbag(fNumData); // Initialization includes filling with 0 // bootstrap aggregating (bagging) -> sampling with replacement: @@ -221,14 +215,7 @@ // {0,1,...,fNumData-1}, which is the set of the index numbers of // all events in the training sample - - for (Int_t n=0;nRndm()); + for (Int_t n=0; nRndm()*fNumData); classpopw[fHadTrue[k]]+=fWeight[k]; @@ -241,7 +228,11 @@ // In bagging procedure ca. 2/3 of all elements in the original // training sample are used to build the in-bag data - datsortinbag=fDataSort; - - ModifyDataSort(datsortinbag,ninbag,jinbag); + TArrayI datsortinbag=fDataSort; + Int_t ninbag=0; + + ModifyDataSort(datsortinbag, ninbag, jinbag); + + const TMatrix &hadrons=fHadrons->GetM(); + const TMatrix &gammas =fGammas->GetM(); // growing single tree @@ -258,39 +249,38 @@ // determined from oob-data is underestimated, but can still be taken as upper limit. - TVector event(fNumDim); - for(Int_t ievt=0;ievt0)continue; - event=TMatrixRow(hadrons,ievt); - fHadEst[ievt]+=fRanTree->TreeHad(event); + for (Int_t ievt=0;ievt0) + continue; + fHadEst[ievt] += fRanTree->TreeHad(hadrons, ievt); fNTimesOutBag[ievt]++; } - for(Int_t ievt=0;ievt0)continue; - event=TMatrixRow(gammas,ievt); - fHadEst[fNumHad+ievt]+=fRanTree->TreeHad(event); + for (Int_t ievt=0;ievt0) + continue; + fHadEst[fNumHad+ievt] += fRanTree->TreeHad(gammas, ievt); fNTimesOutBag[fNumHad+ievt]++; } Int_t n=0; - fErr=0; - for(Int_t ievt=0;ievtGetM(); - TMatrix gammas=fGammas->GetM(); + const TMatrix &hadrons=fHadrons->GetM(); + const TMatrix &gammas=fGammas->GetM(); for (Int_t j=0;jdue to bagging) fDataSort void CreateDataSort(); - void ModifyDataSort(TArrayI &datsortinbag,Int_t ninbag,TArrayI &jinbag); + void ModifyDataSort(TArrayI &datsortinbag, Int_t ninbag, const TArrayI &jinbag); public: Index: trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc (revision 2295) +++ trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc (revision 2296) @@ -84,5 +84,5 @@ // Needs: // - MatrixGammas [MHMatrix] -// - MatrixHadrons {MHMatrix] +// - MatrixHadrons [MHMatrix] // - MHadronness // - all data containers used to build the matrixes Index: trunk/MagicSoft/Mars/mranforest/MRanTree.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanTree.cc (revision 2295) +++ trunk/MagicSoft/Mars/mranforest/MRanTree.cc (revision 2296) @@ -76,5 +76,5 @@ } -void MRanTree::GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue, +void MRanTree::GrowTree(const TMatrix &mhad, const TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue, TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag, TArrayF &winbag,TArrayF &weight) @@ -88,10 +88,4 @@ for (Int_t n=0;nRndm()); + Int_t mvar=Int_t(gRandom->Rndm()*mdim); // Gini index = rrn/rrd+rln/rld @@ -197,12 +194,9 @@ rln=0; rld=0; - wl.Reset(); - - for (Int_t j=0;j<2;j++) - { - wr[j]=tclasspop[j]; - } - - critvar=-1.0e20; + + TArrayF wl(2); // left node + wr = tclasspop; + + Double_t critvar=-1.0e20; for(Int_t nsp=ndstart;nsp<=ndend-1;nsp++) @@ -223,5 +217,5 @@ if (datarang[mvar*numdata+nc]1.0e-5) + if (TMath::Min(rrd,rld)>1.0e-5) { crit=(rln/rld)+(rrn/rrd); @@ -309,6 +303,5 @@ TArrayI &nodepop,TArrayI &nodestart,TArrayF &tclasspop, Int_t nrnodes,TArrayI &idmove,TArrayI &ncase,TArrayI &parent, - TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr,TArrayF &wc, - TArrayF &wl,Int_t ninbag) + TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t ninbag) { // Buildtree consists of repeated calls to two void functions, FindBestSplit and MoveData. @@ -326,18 +319,11 @@ // returns. - Int_t msplit,nbest,ndendl,nc,jstat,ndend,ndstart; + Int_t msplit, nbest, ndendl; Float_t decsplit=0; - Float_t popt1,popt2,pp; - TArrayF classpop; - TArrayI nodestatus; - - nodestatus.Set(nrnodes); - classpop.Set(2*nrnodes); - - nodestatus.Reset(); + TArrayF classpop(2*nrnodes); + TArrayI nodestatus(nrnodes); + nodestart.Reset(); nodepop.Reset(); - classpop.Reset(); - for (Int_t j=0;j<2;j++) @@ -357,13 +343,12 @@ // initialize for next call to FindBestSplit - ndstart=nodestart[kbuild]; - ndend=ndstart+nodepop[kbuild]-1; + const Int_t ndstart=nodestart[kbuild]; + const Int_t ndend=ndstart+nodepop[kbuild]-1; for (Int_t j=0;j<2;j++) tclasspop[j]=classpop[j*nrnodes+kbuild]; - jstat=FindBestSplit(datasort,datarang,hadtrue,mdim,numdata, - ndstart,ndend,tclasspop,msplit,decsplit, - nbest,ncase,jinbag,iv,winbag,wr,wc,wl, - kbuild); + const Int_t jstat=FindBestSplit(datasort,datarang,hadtrue,mdim,numdata, + ndstart,ndend,tclasspop,msplit,decsplit, + nbest,ncase,jinbag,iv,winbag,kbuild); if(jstat==1) { @@ -392,6 +377,6 @@ for (Int_t n=ndstart;n<=ndendl;n++) { - nc=ncase[n]; - Int_t j=hadtrue[nc]; + const Int_t nc=ncase[n]; + const Int_t j=hadtrue[nc]; classpop[j*nrnodes+ncur+1]+=winbag[nc]; } @@ -399,6 +384,6 @@ for (Int_t n=ndendl+1;n<=ndend;n++) { - nc=ncase[n]; - Int_t j=hadtrue[nc]; + const Int_t nc=ncase[n]; + const Int_t j=hadtrue[nc]; classpop[j*nrnodes+ncur+2]+=winbag[nc]; } @@ -410,16 +395,17 @@ if (nodepop[ncur+1]<=fNdSize) nodestatus[ncur+1]=-1; if (nodepop[ncur+2]<=fNdSize) nodestatus[ncur+2]=-1; - popt1=0; - popt2=0; + + Double_t popt1=0; + Double_t popt2=0; for (Int_t j=0;j<2;j++) { - popt1+=classpop[j*nrnodes+ncur+1]; - popt2+=classpop[j*nrnodes+ncur+2]; + popt1+=classpop[j*nrnodes+ncur+1]; + popt2+=classpop[j*nrnodes+ncur+2]; } for (Int_t j=0;j<2;j++) { - if (classpop[j*nrnodes+ncur+1]==popt1) nodestatus[ncur+1]=-1; - if (classpop[j*nrnodes+ncur+2]==popt2) nodestatus[ncur+2]=-1; + if (classpop[j*nrnodes+ncur+1]==popt1) nodestatus[ncur+1]=-1; + if (classpop[j*nrnodes+ncur+2]==popt2) nodestatus[ncur+2]=-1; } @@ -446,5 +432,5 @@ { fNumEndNodes++; - pp=0; + Double_t pp=0; for (Int_t j=0;j<2;j++) { @@ -482,9 +468,35 @@ const Int_t m=fBestVar[kt]; - kt = event(m)<=fBestSplit[kt] ? fTreeMap1[kt] : fTreeMap2[kt]; } return fBestSplit[kt]; +} + +Double_t MRanTree::TreeHad(const TMatrixRow &event) +{ + Int_t kt=0; + // to optimize on storage space node status and node class + // are coded into fBestVar: + // status of node kt = TMath::Sign(1,fBestVar[kt]) + // class of node kt = fBestVar[kt]+2 (class defined by larger + // node population, actually not used) + // hadronness assigned to node kt = fBestSplit[kt] + + for (Int_t k=0;k